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AN  ANALYSIS  OF  RECOVERABLE  ITEM  INVENTORY  SYSTEMS 
WITH  SERVICE  FACILITIES  SUBJECT  TO  BREAKDOWN 

Peter  L.  Knepell,  Ph.D. 

Cornell  University  1979 

The  purpose  of  this  study  Is  to  analyze  an  inventory /maintenance 
system  for  recoverable  Items,  that  Is,  Items  which  are  subject  to 
repair  when  they  fail.  The  repair  of  items  is  performed  by  a main- 
tenance facility  which  has  a fixed  number  of  service  stations  or  chan- 
nels which  are  also  subject  to  failure.  When  an  item  fails,  a demand  is 
immediately  placed  for  a like  replacement  from  a spare  pool.  The  failed 
part  is  sent  to  the  repair  facility  to  be  serviced  on  a first-come, 
first-served  basis.  The  spare  pool  is  replenished  when  repair  on  the 
item  is  completed.  When  a service  station  fails,  repair  is  initiated 
immediately  and  the  failed  server  is  replaced  by  an  operative  spare 
server  if  one  is  available.  This  analysis  is  limited  to  a single- 
echelon system  with  no  outside  sources  of  supply  or  repair. 

The  objective  of  this  study  is  to  model  the  system  described  in 
order  to  observe  the  relationship  of  system  performace  to  spare  stock 
levels  and  service  facility  design.  Specifically,  the  model  is  used  to 
minimize  the  total  expected  unit  backorders  given  an  investment  con- 
straint on  the  number  of  spare  items,  service  channels  and  spare  servers 
in  the  system.  For  long  range  planning  purposes,  this  is  accomplished 
for  a system  with  demands  which  are  stochastic  and  stationary  in  nature. 
An  extension  is  provided,  to  consider  the  case  where  demands  are  non- 
statlonary  and/or  the  time  dependent  behavior  of  the  system  needs  to  be 


described 


In  order  to  express  the  total  expected  unit  backorders,  a repre- 
sentation for  the  distribution  of  the  number  of  units  requiring  repair 
is  needed.  Approximations  are  developed  using  diffusion  techniques 


since  the  actual  distributions  are  difficult  to  express.  The  diffusion 
approximation  is  applied  to  an  optimization  problem  to  provide  the  best 
allocation  of  investments  in  the  system.  A simple  solution  algorithm  is 
given. 

Finally,  a view  of  the  time-dependent  behavior  of  the  system  is 
provided.  The  problem  is  decomposed  into  finding  the  distributions  for 
(1)  the  number  of  units  in  requiring  repair  given  no  service  channel 
failures  and  (2)  the  time  between  service  channel  failures.  We  provide 
a brief  review  of  the  literature  for  the  first  distribution  and  an  in- 
depth  study  of  the  latter  distribution.  ° 
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CHAPTER  I 


INTRODUCTION 


The  purpose  of  this  study  is  to  analyze  an  inventory /maintenance 
system  for  recoverable  items,  that  is,  items  which  are  subject  to 
repair  when  they  fail.  The  repair  of  items  is  performed  by  a main- 
tenance facility  which  has  a fixed  number  of  service  stations  or  chan- 
nels which  are  also  subject  to  failure.  When  an  item  fails,  a demand  is 
immediately  placed  for  a like  replacement  from  a spare  rcol.  The  failed 
part  is  sent  to  the  repair  facility  to  be  serviced  on  a first-come, 
first-served  basis.  The  spare  pool  is  replenished  when  repair  on  the 
item  is  completed.  When  a service  station  fails,  repair  is  initiated 
immediately  and  the  failed  server  is  replaced  by  an  operative  spare 
server  if  one  is  available.  This  analysis  is  limited  to  a single- 
echelon system  with  no  outside  sources  of  supply  or  repair. 

The  objective  of  this  study  is  to  model  the  system  described  in 
order  to  observe  the  relationship  of  system  performace  to  spare  stock 
levels  and  service  facility  design.  Specifically,  the  model  is  used  to 
minimize  the  total  expected  unit  backorders  given  an  investment  con- 
straint on  the  number  of  spare  items,  service  channels  and  spare  servers 
in  the  system.  For  long  range  planning  purposes,  this  is  accomplished 
for  a system  with  demands  which  are  stochastic  and  stationary  in  nature. 
An  extension  is  provided,  to  consider  the  case  where  demands  are  non- 
stationary and/or  the  time  dependent  behavior  of  the  system  needs  to  be 
described. 

These  inventory /maintenance  systems  involve  vas  : capital  invest- 
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raents;  hence,  the  design  and  control  of  these  systems  are  a great 
concern  for  managers.  In  large  scale  industrial  and  military  activi- 
ties, a majority  of  the  inventory  items  are  inexpensive  consumable  (non- 
recoverable)  units;  however,  a large  proportion  of  the  inventory  invest- 
ment is  for  spare  stock  levels  of  recoverable  items.  Sherbrooke  [32] 
states  that  recoverable  item  spares  in  the  Air  Force  account  for  78 
percent  of  the  total  investment,  amounting  to  approximately  five  billion 
dollars.  Currently,  automated  repair  stations,  costing  up  to  16  million 
dollars  each,  are  being  purchased  by  military  organizations  to  repair 
units  which  cost  an  average  of  100  thousand  dollars  each. 

To  provide  a better  understanding  of  the  structure  of  this  system, 
we  will  describe  a specific  example  where  the  units  which  fail  are 
sophisticated  aircraft  electronic  (avionics)  components,  such  as  radar, 
navigation  instruments  and  radios.  Spare  units  are  stocked  at  the 
airfield  where  the  demands  occur.  The  unit  failure  rates  are  usually 
low  and  failures  occur  independently.  Occasionally,  a rash  of  break- 
downs occur  in  a particular  item  or  the  failure  of  one  item  may  induce 
the  failure  of  different  units.  These  sources  of  dependent  demand  are 
infrequent  and  difficult  to  predict  and,  therefore,  are  not  considered 
in  the  analysis.  In  practice,  a small  number  of  units  are  sent  to 
another  facility  or  higher  echelon  level  for  repair  or  replacement. 

The  repair  stations  are  situated  at  the  airfield.  They  also  in- 
volve sophisticated  electronic  components  and  their  failure  character- 
istics are  similar  to  those  of  the  recoverable  items.  In  addition, 
these  stations  are  periodically  out  of  service  for  modification,  pre- 
ventive maintenance,  and  calibration.  If  a service  station  fails,  a 
high  priority  is  placed  on  its  repair  since  service  interruptions 
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ultimately  affect  the  number  of  operational  aircraft.  A supply  of  spare 
components  for  the  repair  stations  is  usually  provided. 

Thus  the  objective  in  designing  the  system  in  this  example  is  to 
provide  the  optimal  investment  allocation  for  spare  units  and  service 
facilities  at  the  airfield  so  that  the  maximum  number  of  aircraft  are 
operationally  ready. 

We  begin  the  study  in  Chapter  II  with  a brief  description  of  some 
previously  developed  models  of  recoverable  item  inventory  systems.  This 
is  followed  by  a justification  for  the  use  of  the  expected  number  of 
unit  backorders  as  a measure  of  system  performance.  After  some  basic 
assumptions  are  listed,  a mathematical  statement  of  our  model  is  pro- 
vided. 

Chapters  III  and  IV  provide  a planning  tool  for  managers  to  use 
when  designing  a recoverable  item  inventory  system.  Chapter  III  devel- 
ops methods  for  obtaining  the  stationary  probability  distribution  of  the 
number  of  units  being  repaired.  We  use  this  distribution  to  compute  the 
expected  number  of  unit  backorders.  Approximations  are  developed  using 
diffusion  techniques  since  the  actual  distributions  are  difficult  to 
express.  In  Chapter  IV,  the  diffusion  approximation  is  applied  to  an 
optimization  problem  to  decide  on  the  allocation  of  investments  in  the 
system.  A simple  solution  algorithm  is  given. 

A view  of  the  time  dependent  behavior  of  the  system  is  given  in 
Chapter  V.  The  problem  is  decomposed  into  finding  the  distributions  for 
(1)  the  number  of  units  requiring  repair  given  no  service  channel 
failures  and  (2)  the  time  between  service  channel  failures.  We  provide 


a brief  review  of  the  literature  for  the  first  distribution  and  an  in- 
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depth  study  of  the  latter  distribution.  The  final  chapter  contains  some 
closing  remarks  and  suggestions  for  future  research. 


CHAPTER  II 


THE  MODEL 

Much  attention  has  been  focused  on  attempts  to  model  inventory 
systems  like  the  one  described  earlier.  Feeney  and  Sherbrooke  [5] 
examined  single-echelon  recoverable  item  systems  where  demand  was 
generated  by  a compound  Poisson  process.  Sherbrooke  [32]  extended  the 
results  to  a two  echelon  system  in  a model  he  called  METRIC  (Multi- 
Echelon  Technique  for  Recoverable  Item  Control).  Muckstadt  [23]  ex- 
tended the  METRIC  model  to  include  part  hierarchies.  All  of  these 
papers  assumed  that  the  service  facility  has  adequate  capacity  to  repair 
all  items  without  delay  (i.e.,  the  infinite  server  assumption).  Gross, 
et.  al.,  [2,9]  considered  a recoverable  item  system  with  finite  service 
capacity.  They  modeled  their  system  as  a classical  machine-repairman 
problem.  However,  they  too  assumed  the  servers  were  reliable. 

Typically,  service  facilities  are  constrained  in  their  capacity  to 
a finite  number  of  servers  and,  in  some  cases,  the  servers  are  subject 
to  failure.  Under  these  considerations  it  is  important  to  consider  the 
effects  of  service  facility  design  on  overall  system  performance.  Some 
of  the  work  done  in  this  area  will  be  discussed  later;  in  general,  very 
few  results  are  available  in  the  context  of  production-inventory  con- 
trol, 


2.1  Purpose  of  the  Model 

A model  for  a single-echelon,  recoverable  item  inventory  system 
will  be  developed  in  this  chapter.  The  model  can  then  be  used  to 
quantify  the  relationships  between  CD  service  reliability  and  capacity 
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and  (2)  overall  system  performance.  While  it  will  be  useful  as  a design 
tool  for  managers,  it  is  not  intended  to  help  them  make  day-to-day 
decisions  in  the  dynamic  environment  of  the  inventory  system.  Although 
it  is  our  objective  to  create  a realistic  model,  some  simplifying 
assumptions  will  be  made  to  facilitate  the  analysis.  The  most  important 
step  is  to  establish  a meaningful  performance  measure. 


2.2  The  Expected  Backorder  Objective  Function 

We  shall  use  the  sum  of  the  expected  unit  backorders  as  our  per- 
formance measure.  Consider  the  system's  operation  for  a fixed  number  of 
days,  and  count  the  total  number  of  days  in  which  units  are  backordered 
for  that  period.  The  expected  value  of  this  number  divided  by  the 
number  of  days  in  the  period  gives  the  expected  backorders  per  day.  Our 
goal  is  to  minimize  this  function.  Note  that  by  this  definition,  a ten 
day  backorder  is  equivalent  to  ten  backorders  for  one  day. 

Other  performance  measures,  such  as  NORS  rate,  fill  rate,  and  ready 
rate,  are  not  as  versatile  as  expected  backorders  when  modeling  recover- 
able item  inventory  systems.  In  Air  Force  parlance,  the  NORS  (not 
operationally  ready,  supply)  rate  is  considered  an  excellent  measure  of 
logistics  support.  This  figure  represents  the  minimum  number  of  air- 
craft which  cannot  perform  a mission  due  to  supply  backorders.  It  has 
the  advantage  of  measuring  the  direct  impact  that  the  inventory  system 
has  on  the  fleet  it  supports.  Unfortunately,  it  is  a difficult  measure 
to  quantify  and  use  in  inventory  models.  For  example,  if  ten  different 
aircraft  are  grounded,  each  due  to  a different  component  being  back- 
ordered, then  only  one  aircraft  is  considered  NORS.  This  is  further 


complicated  by  component  interchangeability,  substitutability  and 
redundancy.  Once  quantified,  the  NORS  function  is  not  separable  and, 
hence,  is  difficult  to  work  with. 

Another  measure,  fill  rate,  is  defined  as  the  fraction  of  demands 
that  are  immediately  satisfied  by  supply.  This  measure  ignores  the 
length  of  time  a unit  is  backordered.  Sherbrooke  [32]  points  out  that 
when  fill  rate  is  employed  in  a multi-echelon  inventory  system,  managers 
are  encouraged  to  concentrate  nearly  all  stock  at  the  lowest  echelon. 
While  backorders  will  be  infrequent,  they  will  have  long  durations. 
Another  disadvantage  is  that  a "fill"  is  usually  defined  as  an  immediate 
satisfaction  of  a demand.  If  a short  delay  in  satisfaction  is  accept- 
able, the  resulting  optimal  policy  may  be  considerably  different.  In 
fact,  as  longer  delays  are  accepted,  the  more  closely  the  results 
resemble  those  of  the  backorder  criterion. 

Ready  rate  is  defined  as  the  fraction  of  items  which  are  not  back- 
ordered. This  measure  does  not  reflect  the  number  of  units  backordered 
on  a particular  item.  Thus,  it  is  conceivable  that  inexpensive  items 
will  be  stocked  heavily  in  favor  of  the  expensive  items.  Then  back- 
orders will  accumulate  only  on  the  relatively  few  expensive  items  and 
the  system  performance,  measured  by  ready  rate,  will  appear  excellent 
while  large  numbers  of  backorders  exist  for  expensive  items. 

The  expected  backorder  criterion  combines  the  number  of  backorders 
and  the  length  of  each  backorder  as  a penalty.  It  eliminates  the  need 
to  determine  arbitrary  backorder  and  holding  costs  and  provides  a direct 
measure  of  support  that  the  inventory  system  provides.  Sherbrooke  [32] 
mentions  a single-echelon  example  in  which  fill  rate,  ready  rate,  and 
expected  backorder  objective  functions  provide  essentially  identical 
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stockage  policies.  However,  when  applied  to  multi-echelon  problems,  the 
expected  backorder  criterion  yields  more  reasonable  results.  Addition- 
ally, the  expected  backorder  function  is  convex  and  separable,  proper- 
ties which  are  computationally  helpful  and  not  necessarily  possessed  by 
other  criteria. 

2.3  Basic  Assumptions 

A list  of  the  basic  assumptions  made  for  this  model  will  be  given, 
followed  by  an  expanded  discussion  of  each. 

1.  The  demand  process  for  each  of  m different  items  is  a Poisson 
process.  All  demands  occur  independently  at  a rate  i=l,,..,m. 

2.  With  each  demand,  units  are  exchanged  on  a one-for-one  basis. 

3.  All  units  turned  in  are  serviced. 

4.  Service  times  are  stochastically  independent  and  exponentially 
distributed  at  rate  y.  There  are  at  most  C service  channels  available. 
If  service  is  interrupted  by  a channel  failure,  then  the  unit  is  immedi- 
ately moved  to  the  next  available  service  channel  and  service  is  resumed 
without  delay. 

5.  There  is  no  batching  of  items  for  repair.  Items  are  serviced 
on  a first-in,  first-served  basis. 

6.  Service  channels  fail  independently  as  Poisson  events  at  a 
rate  £. 

7.  Pailed  servers  are  repaired  immediately  with  exponentially 
distributed  repair  times  at  a rate  r). 

8.  Failed  channels  are  replaced  Instantaneously,  if  a spare  is 
available.  If  no  spare  is  available,  channels  are  replaced  in  order  of 
breakdown  times  when  repaired  servers  become  available.  L spare  servers 
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are  provided. 

The  first  assumption  implies  that  the  arrival  rate  of  units  does 
not  depend  on  the  size  of  the  population.  In  a standard  application, 
there  is  a finite  source  of  demand;  however,  when  looking  at  the  scenarios 
we  are  modeling,  this  assumption  is  valid.  For  example,  consider  a 
fleet  of  aircraft  as  generating  units  requiring  service.  As  aircraft 
units  are  backordered,  the  number  of  operational  aircraft  decreases. 

Since  the  flying  schedule  is  fixed,  the  remaining  aircraft  will  have  to 
fly  more  to  satisfy  the  schedule.  Since  we  assume  an  aircraft's  failure 
rate  is  directly  related  to  its  usage,  the  perceived  fleet  failure  rate 
is  assumed  to  remain  the  same.  This,  of  course,  neglects  the  possibil- 
ity that  a large  proportion  or,  for  that  matter,  the  entire  fleet  could 
be  grounded  at  the  same  time.  In  a later  chapter  we  will  allow  the 
demand  rate  to  vary  over  time  to  reflect  changing  flying  schedules  or 
failure  characteristics. 

The  second  assumption  reflects  an  (S,S-1)  inventory  policy.  The 
next  assumption  states  that  the  system  is  conservative  (i.e.,  no  condem- 
nations). In  practice,  the  inventory  items  are  expensive,  so  these 
assumptions  reflect  a reasonable  policy. 

The  fourth  assumption  gives  non-preemptive  priority  to  a unit  whose 
service  is  interrupted.  Theoretically,  given  the  memoryless  property  of 
the  exponential  distribution,  this  assumption  does  not  matter.  It  will 
become  apparent  later  that  exponentially  distributed  interarrival  times 
are  not  necessary  for  the  approximation  techniques  U9ed.  This  assump- 
tion is  made  because  the  approximation  method  proposed  was  tested 
against  systems  with  exponentially  distributed  service  times. 

The  sixth  assumption  implies  that  service  channels  can  fail  even 


when  idle.  This  is  realistic  since  calibration  tests,  preventive  main- 
tenance and  modifications  are  typical  for  the  systems  being  modeled. 

The  next  assumption  implies  that  there  is  an  adequate  number  of  repair- 
men to  work  on  the  failed  channels.  This  can  be  altered  to  specify  a 
limited  number  of  repairmen;  however,  this  adds  to  the  notational  and 
computational  complexity,  but  does  not  alter  the  method  of  analysis. 

In  Chapter  I,  it  was  mentioned  that  the  service  channels  are  large 
and  costly  to  establish.  Spares  for  these  channels  are  relatively 
inexpensive.  Therefore,  it  is  reasonable  to  assume  that  while  C + L 
servers  are  provided,  only  C are  useable  service  channels  and  the 
remaining  L must  be  held  in  reserve. 

2.4  Mathematical  Statement  of  the  Model 

The  objective  of  inventory  managers  is  to  provide  the  greatest 
system  performance  given  a fixed  budget.  Thus,  our  goal  is  to  minimize 
total  expected  backorders  outstanding  at  any  point  in  time  subject  to  a 
budget  constraint.  This  model  will  be  different  from  those  mentioned 
earlier  because  the  Investment  constraint  links  the  purchase  of  unit 
spares,  service  channels,  and  spare  servers.  Thus  the  system  perform- 
ance will  be  dependent  on  the  service  facility  design  as  well  as  the 
allocation  of  funds  to  unit  spares. 

Suppose  p(n|*)  represents  the  probability  that  n units  are  in  the 
service  facility  (in  service  or  awaiting  service).  We  know  this  number 
will  be  a function  of  the  input  and  the  output  of  the  service  facility. 
The  parameter  characterizes  the  input  process  for  item  i and  the 
parameters  y C,  and  L determine  the  output  process.  For  each  unit 
of  type  i we  can  express  the  expected  number  of  backorders  outstanding 
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at  any  time  as 

l (x  - s ) p(x|A  ,y,S,n,C,L) . (2.1) 

x>s1 

Backorders  in  some  items  may  be  considered  more  serious  than  for  others. 
In  this  case,  we  can  weight  the  backorder  function  in  (2.1)  with  an 
essentiality  factor,  E^. 

The  costs  of  unit  spares  and  service  channels  are  considered  to  be 
linear.  Since  the  initial  setup  cost  for  service  channels  is  very 
large,  there  will  be  different  costs  for  service  channels  and  their 
spares.  An  expression  for  the  investment  constraint  is 

m 

C*C  + L*C  + y C.s.  < B.  (2.2) 

c s ^ i i - 

If  we  are  concerned  with  a long  range  planning  tool,  we  must  be 
careful  to  design  a service  facility  which  can  provide  adequate  service 
over  an  extended  period  of  time.  Clearly,  the  potential  output  rate 

must  be  at  least  as  large  as  the  input  rate: 

m 

A = l A < yC,  (2.3) 

i=l 

where  C is  the  expected  number  of  service  channels  operational  at  any 
point  in  time.  The  mathematical  derivation  of  this  constraint  will  be 
provided  later. 

Combining  all  the  statements  above  we  have  a mathematical  statement 
of  the  model  as  follows: 

m 

minimize  £ £ (x  - s^)  p(x| A^,y ,s*n,C,L) 

i-1  x>sL 

m 

CC  + LC  + 7 C.s.  < B, 

c s i i- 


subject  to 


(P) 
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and  X < yC, 

where  C,  L,  and  Sj^  are  non-negacive  integers,  1-1,... m. 

For  future  reference,  this  optimization  problem  will  be  denoted  as 
problem  P. 


CHAPTER  III 

STATIONARY  DISTRIBUTION  ANALYSIS 

To  apply  the  model  developed  in  Chapter  II  as  a long  range  planning 
tool,  we  need  to  find  an  accurate  expression  for  p(n|*),  the  stationary 
distribution  for  the  number  of  units  in  the  system.  To  do  this,  we  view 
the  single-echelon,  recoverable  item  system  as  a queueing  system  in 
which  the  service  facility  has  a finite  number  of  servers,  each  subject 
to  failure.  The  servers  and  their  repair  facility  will  be  another 
queueing  system  imbedded  in  the  first  one.  This  structure  will  be 
exploited  to  develop  the  stationary  probability  distributions  of  units 
in  the  service  system.  Some  analytic  results  will  be  given  and  a 
general  approximation  method  will  be  derived. 

3.1  The  Queueing  Model 

This  section  will  develop  a queueing  model  to  represent  the  system 
introduced  in  Chapter  I.  It  will  become  apparent  that  this  system  has 
much  in  common  with  the  classical  machine-repairman  problem.  The  assump- 
tions given  previously  establish  the  queue  service  discipline  and  allow 
for  the  creation  of  a state  space  for  a continuous-time  Markov  process. 

The  system  to  be  modeled  can  be  described  schematically  (see  Figure 
3.1).  Units  demanding  service  come  from  an  infinite  source.  They  enter 
a multi-server  queue  and  are  served  on  a first-come,  first-served  basis. 
As  soon  as  a unit  joins  the  queue,  a replacement  unit,  if  one  is  avail- 
able, is  immediately  returned  to  the  source  from  the  spare  pool.  If  no 
replacement  is  available  it  is  backordered. 

13 
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Figure  3.1:  Queueing  System  Model 

Servers  are  subject  to  failure  and  are  replaced  by  spare  servers,  if  one 
is  available.  Figure  3.1  illustrates  why  this  system  could  be  called  "a 
two-dimensional  machine-repairman  system." 

We  must  assume  this  system  is  non-saturated  in  order  to  guar  '.ee 
the  existence  of  a stationary  probability  distribution.  To  assure  this, 
as  mentioned  earlier,  we  will  require 

X < HC,  (3.1) 
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then  from  Palm's  Theorem  we  know  that  pn  is  Poisson  with  rate  X/u  . If 
Ol  and  there  are  no  service  station  failures  (i.e.,  £>0  or  L-**>  ),  then 
Pn  is  geometric  with  parameter  X/u,  since  we  have  a simple  M/M/1 
queue.  If  there  are  no  spare  service  stations  (L-0)  and  C is  finite, 
then  p^  can  be  determined  using  the  results  of  Mitrany  and  Avi-Itzhak  [20]. 

Since  the  servers  can  fail  at  any  time,  the  number  of  operational 
channels  is  independent  of  the  number  of  units  in  the  system.  (The 
reverse,  however,  is  not  true.)  Thus,  the  stationary  probability 
distribution  of  the  number  of  servers  operational  can  be  obtained 
separately.  This  subsystem  can  be  viewed  as  a classical  machine-repair- 
man problem  where  the  servers  and  their  spares  are  the  "machines"  and 
there  are  always  an  adequate  number  of  "repairmen." 


Theorem  3.1. 

If  we  have  a system  with  C service  channels  subject  to  failure,  L 
spare  servers,  and  C + L repairman  where: 

i)  the  servers  fail  independently,  as  Poisson  events  at  rate  C, 

and 

ii)  the  repair  times  are  independent  and  exponentially  distributed 
with  rate  n, 

then  the  stationary  probability  of  having  k operational  servers,  q^, 
exists  and  is  given  by 


16 


where 


'C+L 


L k 

k °tL  L 

- 

= 

I 

ft 

k + y cLc: 

£)k 

k=0  k! 

k 

k=L+l  (C+L-k)  .'k! 

nj 

-1 


(3.3) 


A simple  proof  of  these  results  can  be  found  in  reference  6. 

Thus  the  stationary  probability  of  having  k repair  channels  oper- 
ational is 


, 0 < k < C - 1, 


Pr{k  channels  up}  =s 


C+L 

l \ 

n=C 


k = C. 


(3.4) 


In  the  special  case  where  L«0  (i.e.,  no  spares)  we  have 

i C-k 


qk  = 


fi]  q 

nj  4 


C+L 


where 


0 £ k <_  C , 

1 


HC+L  C 

I 

j-  k=0 

Defining  P = — we  have 


i + £! 

nj 


(3.5) 

(3.6) 


qk  = 


(1  + P) 


1 + P 


l1  + Pj 


C-k 


, 0 £ k £ C. 

(3.7) 


So  when  no  spare  servers  are  provided,  the  stationary  probability 
distribution  is  binomial. 

The  conditional  probability  that  k servers  are  operational  given 

that  a failure  is  just  about  to  occur  is  not  the  same  as  since  the 

failed  service  units  come  from  a finite  source.  A derivation  for  this 
» 

new  distribution,  q^,  can  be  found  in  reference  9. 
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Corollary  3.1.  Let  the  same  conditions  as  Theorem  3.1  hold.  Then  the 
stationary  probability 

• Pr  { k servers  operational  | a server  failure  Is  about  to  occur) 
exists  and  is  given  by 




C 

C - l (C-n)q 
n-0 


J*k 

C 

C - l (C-n)q 
n-0 


, 0 < k < C 


, C +1  _<  k £ L+C. 

(3.8) 


3.2  Analytic  Results 

The  typical  procedure  used  to  derive  the  line  length  distribution 
involves  probability  generating  functions.  A Markovian  state  space  is 
designed,  balance  equations  are  developed  and  a generating  function  is 
derived.  The  poles  of  the  generating  function  provide  information  which 
is  otherwise  difficult  to  derive.  The  previous  work  found  in  the 
literature  has  been  limited  to  cases  where  no  spare  service  channels 
are  available.  This  work  will  be  discussed  in  greater  detail  and  then  a 
derivation  for  the  "simple"  case  of  one  server  and  one  spare  will  be 
given. 


3.2.1  Previous  Results 

The  earliest  published  work  done  on  queueing  systems  with  service 
station  subject  to  breakdown  was  by  White  and  Christie  in  1958  [34], 

The  most  comprehensive  article  on  single  server  queueing  systems  with 
Poisson  inputs  and  exponential  service  times  was  written  in  1963  by  Avi- 
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Itzhak  and  Naor  [1],  They  allowed  a general  distribution  for  server 
breakdowns  and  repair  and  were  only  able  to  derive  the  expected  line 
length  and  expected  waiting  time.  Shogan  in  1977  provided  the  line 
length  distribution  for  a single  server  system  [33].  His  results  are 
summarized  below. 


Theorem  3.2  (Shogan) 

Suppose  we  have  a system  with: 

(i)  independent,  exponentially  distributed  inter-arrival  times, 
rate  A, 

(ii)  independent  Erlang  distributed  service  times  with  mean  1/y 
and  shape  parameter  k, 

(iii)  a single  server,  no  spares,  with  exponential  inter-failure 
times,  rate  £,  and 

(iv)  Erlang  distributed  repair  times  with  mean  1/n  and  shape 
parameter  m. 


If  A/y  < n/(C 


P - Pr 
riJ 


exists  and  is 
P00 
Pm,  j 

Pi,j 


0,j+l 


+ n) , then  the  stationary  probability  distribution 
{ i phases  of  repair  required  on  the  server, 
j phases  of  service  in  the  system) 


defined  by 

the  recursive 

equations  (for  j ■ 

0,1,2,. 

..): 

- n/(£+n)  - 

■ A/y, 

(3.9) 

= (A-mn)  ^ 

(£?n  , + AP 

0 3 J ® 9 

j~k)  * 

(3.10) 

* (A-mn)  ^ 

(mnPi+i , j + xpi>J-k)-  i— lj  "'2 

* • • • 

(3.11) 

j 

j j 

- (kM)_1  (A  l ?nn 

♦ 5 l P0  - »n  I 

?^  * 

1 f n 

(3.12) 

n-J+l-k  ’ 

n-0  ’ n=0 
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where  Is  zero  If  It  has  a negative  subscript  and  summations  are  zero 
If  the  lover  limit  of  summation  Is  negative. 


The  expressions  In  Theorem  3.2  are  typical.  A closed  from  solution 
for  the  line  length  probabilities  has  not  been  found.  The  recursive 
relations  are  a direct  result  of  the  balance  equations  and  Pqq  was 
obtained  from  the  pole  of  a lengthy  generating  function.  After  solving 
for  all  the  state  probabilities,  the  steady  state  probability  of  having 
n customers  in  the  system  is  given  by: 


m 

I Pi  0 

i-0  ’ 


for  n-0 


m nk 

l l Pii 

i-0  j*(n-l)k+l  J 


for  n > 0. 


(3.13) 


The  groundwork  for  multi-server  queues  with  server  breakdowns  and 
no  spares  has  been  established  in  an  article  by  Mitrany  and  Avi-Itzhak 
in  1968  [20].  Their  results  for  a two  server  system  are  summarized 

below. 


Theorem  3.3  (Mitrany  and  Avi-Itzhak) 

Suppose  we  have  a system  with: 

(I)  Independent,  exponentially  distributed  inter-arrival  times, 
rate  X, 

(II)  Independent,  exponentially  distributed  service  times  with 
mean  1/p, 

(ill)  two  servers,  no  spares,  with  Independent,  exponentially 
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distributed  inter-failure  times,  rate  £ , and 

(iv)  independent,  exponentially  distributed  repair  times,  with 
mean  1/n  . 

Define 

z - [ (A+y+C+n)  - / (A+y+C+n)2  - 2 Ay' ]/2A  , (3.14) 

N » 2yn  - An  - AC  , and  (3.15) 

D = u(C+n)  [2(C+n) (2u+A+2Cz)  + A(l-z) (2y+A) ] . (3.16) 

If  A/p  < 2n  / (C  +H  ),  then  the  steady  state  probability 

Pij  “ Pr  { i operating  servers,  j units  in  the  system  } 
exists  and  is  defined  by  the  recursive  equations  (for  i - 0,1,2): 


P00  = ^2t4U  + (2X  + 4^)ZJ'N  / (A+n)D  , (3.17) 

P10  - (A-H1)  Pqo  / £ , (3.18) 

P20  3 [ (A+2n  )U  + ( 2£n  -Ay  ) z ] • N / D,  (3 . 19) 

[jy+iC+A+(2-i)n]  =■  (3-i>nPi-i,j  + (i+1)spi+l,j  + APt  i_1 

+ (j+l)ypijj+1  f for  j < i,  (3.20) 

[iy+iC+A+(2-i)n]  - (3-OnP^j  + (i+l)CP1+l  j +AP1j_1 

+ 1UPi,j+l  * for  j ^ i,  (3.21) 


^4 


21 


where  is  zero  if  1 > 2 or  any  subscript  is  negative. 

The  steady  state  probability  of  having  n customers  in  the  system  is 
simply 

The  work  required  for  n > 2 is  computationally  and  notationally 
intractable  and  Mltrany  and  Avl-Itzhak  suggest  numerical  methods.  In  a 
related  article  published  in  1973,  Yechlali  was  also  unable  to  derive 
closed  form  solutions  except  in  very  specific  cases.  He  stated  that, 

"In  general,  no  closed  form  relations  are  available  for  the  probabil- 
ities q,  and,  except  for  numerical  results,  no  analytic  comparison  to 
the  elegant  results  of  the  classical  M/M/1  queue  can  be  made."  [35] 

3.2.2  The  Single  Server,  Single  Spare  System 

This  section  will  show  a technique  to  solve  analytically  for  the 
stationary  line  length  distribution  for  a single  server  queue  subject  to 
breakdown,  where  one  spare  server  is  provided.  In  addition  to  the 
assumptions  stated  in  Section  2.3,  we  will  also  assume  that  operable 
spares  do  not  fall  while  held  in  the  spare  server  pool.  This  is  a 
reasonable  restriction  considering  the  application;  however,  it  can  be 
lifted  without  significantly  affecting  the  analysis. 

The  solution  procedure  to  be  used  follows.  A state  space  for  a 
continuous-time  Markov  process  will  be  defined.  The  Markovian  nature 
of  the  state  space  allows  us  to  describe  the  flows  in  the  system  with 
balance  equations.  These  equations  cannot  be  solved  directly  so  proba- 
bility generating  functions  are  derived.  These  will  assist  us  in 


2 

T P.  for  n > 0. 

i n — 
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solving  for  three  of  the  stationary  probabilities.  When  these  proba- 
bilities are  known,  the  entire  distribution  can  be  obtained  using  the 
balance  equations. 

We  start  by  defining  a state  space  for  this  system  as 

{(i,j)  | i - number  of  operable  servers  ■ 0,1,2, 

j ■ number  of  units  in  the  system  ■ 0,1,...}  . 

Since  the  state  transition  times  are  all  independent  and  exponentially 
distributed,  the  process  described  is  Markovian.  Define  ^ as  the 
steady  state  probability  of  being  in  state  (i,j).  The  transition  flows 
are  illustrated  in  Figure  3.2. 


Figure  3.2:  State  Space  Transitions 
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The  balance  equations  for  this  system  are. 

(x+2n)P0j  - xpq fj-1  + 5Pl#j  ’ j=0>1 

(X+c+n)P10  - 2nP00  + SP20  + vipu  . 

(X+5+u+n)P1;)  - 2nP0j  + ^P2j  + ^pi,j+i  + xpi,j-i  , J-i.2 

(x+C)P20  - hp10  + up21  , 

(x+^u)P2j  - nPLj  + up2J+1  + ’ J-1’3 


Define  the  probability  generating  functions  (pgf) 

G(z)-I  P±j*j  . i“0 » 1 » 2 • 

1 j-0  J 

Note  that  these  generating  functions  evaluated  at  z ■ 1 
stationary  probability  of  the  number  of  operational  servers. 
Theorem  3.1  we  have 


GjU) 


i-0 

9 2 

C + 2nS  + 2n 

» 

2n£ 

. i“l 

S2  + 2nC  + 2nz 

» 

(3.22) 

(3.23) 

(3.24) 

• • • * 

(3.25) 

,...  (3.26) 

give  the 
Applying 
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Also  noee  chat 

G0(l)  + G^l)  + G2C1)  -1. 

We  can  solve  for  each  pgf  by  multiplying  equations  (3.22)  through 
(3.26)  by  the  appropriate  z^  and  summing.  The  result  is: 


[\(l-z)+2n]G0(z)  - CG1(z)  - 0,  (3.29) 

-2nzGQ(z)  + tXz(l-z)+v(z-l)Hz+nz]G1(z)  - 5zG2(z)  - U(z-l)P1(), 

(3.30) 

-HzG1(z)  + [\z(l-z)+u(z-l)+£z]G.,(z)  = u(z-l)p2o. 

(3.31) 

After  some  elementary  row  operations  we  have  an  equivalent  matrix  form 
A(z)'£(z)  - u(z-l)  b,  where: 


A(z)  * 


&(z) 


If  we  had  values  for  pqq»  P^q  and  P2q  we  could  solve  for  all  P^  ^ via 
the  balance  equations.  Equation  (3.22)  provides  one  equation  in  Pqq  and 
P^q.  The  above  pgf’s  will  be  used  to  provide  two  more  independent  equations. 
Using  Cramer's  method  for  solving  simultaneous  equations  we  have: 

| A. (z)  | 

G.  (z)  * — H(z-l)  , i-0,1,2,  (3.32) 

1 |A(z) | 

where  | A^(z)|  is  the  determinate  of  A(z)  with  column  (i  + 1)  replaced 
by  b^.  Examining  row  2 of  A(z)  it  is  clear  that  | A(z)  | ■ (z-1)  J Q(z)  | 
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where 


X(l-z)+2n 


-Xz+u 


-Xz+u 


Xz(l-z)-y (l-z)+£ 


So  we  have 


Gi(z) 


|A±<g) | 

Iq(z)  I 


. 1-0, 1,2. 


(3.33) 


By  its  definition,  G^(z)  must  be  continuous  and  bounded  on  the  interval 

[0,1]  and  we  know  G^l)  - where  q1  is  given  in  (3.27).  Thus  we  have 
G (1) -10(1)1 

»|Ai(l)|  , i « 0,1,2.  (3.34) 

It  is  easy  to  show  that  these  three  equations  are  dependent;  therefore, 

we  only  get  one  useful  equation  in  and  P For  i - 0 we  get 


u(£2+2£n+2n2) 


P + P 
10  20 


(3.35) 


If  we  had  a root  for|Q(z)|on  the  interval  (0,1),  then  by  continuity  of 
Gq(z)  on  (0,1)  we  would  have  another  independent  equation  for  P,  n and 
C‘  Evaluating  |Q(z)|  at  z * 0 and  z = 1 we  have: 


I Q (o) | - -(X  + 2n)  u2  < o 


(3.36) 


lQ(i) I - -2n(X  - U)  (C+  n)  - x^2 

Forcing  I Q(l) | > 0 yields  the  relationship: 


(3.37) 


X < 2Cn  + 2n 
U A , 


S2  + 2£n  + 2n2 


qi  + V 


(3.38) 


p 
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This  is  the  familiar  constraint  requiring  the  expected  service  demanded 
per  unit  time  to  be  less  than  the  expected  service  available.  Given 
that  this  condition  is  met,  the  fourth  degree  polynomial  |Q(z)  | has  a 
root,  z^,  on  the  interval  (0,1).  Then  we  have  |Aq(z^)  | - 0,  or 

equivalently 

[(Xz1-y)(l-z1)+Cz1]P10  + 5zjP20  = 0 . (3.39) 


Equations  (3.22),  (3.35),  and  (3.39)  provide  three  independent  equations 
in  Pqq,  P^O’  and  P20*  Usin8  these,  the  balance  equations  and  the 
traffic  intensity  condition  (3.38)  we  have: 


Theorem  3.4 

Suppose  we  have  a system  with: 

(i)  independent,  exponentially  distributed  inter-arrival  times, 
rate  X, 

(ii)  independent,  exponentially  distributed  service  times  with 
mean  1/y  , 

(iii)  one  server  and  one  spare,  with  independent,  exponentially 
distributed  channel  inter-failure  times,  rate  £ , and 

(iv)  independent,  exponentially  distributed  repair  times  with  mean 


l/n  • 

Define 


K(z) 

D(z) 

qo 


(Xz-y)(l-z)  + £z  , 

[X(z-l)-2n] [Xz-y] [K(z)+nz]  - X£zK(z)  , 

i! 

£2  + 2riS  + 2n2 


(3.40) 

(3.41) 


and 


(3.42) 


26 


If  X/u  < 1 - qg,  then 

a)  D(z)  has  a root,  z^,  on  the  interval  (0,1),  and 

b)  the  stationary  probability  distribution 

■ Pr  ( i servers  are  operative, 

j units  are  in  the  system  } 
exists  and  is  defined  by: 

-q0D(l)  K(z^) 

sVf^  - K(zx)3 

ViD(1) 

» 

[5zl  - K(z^] 

CP10  / (X+2n) 
and  the  balance  equations  (3.22)  through  (3.26). 


20 


10 


00 


(3.43) 


(3.44) 


(3.45) 


3.3  Diffusion  Approximations  For  Queueing  Systems 

It  should  be  evident  that  analytic  results  for  large  systems  are  al- 
gebraically cumbersome  and  do  not  provide  any  insight  into  the  nature  of  the 
desired  distribution.  A computationally  efficient  and  relatively  simple 

i 

approximation  technique  would  greatly  assist  us  in  solving  the  optimization 

I 

problem  proposed  earlier.  For  the  situation  that  we  art  modeling,  diffusion 
approximations  provide  simple  and  accurate  representations  for  the  queue 
size  distribution. 

The  application  of  diffusion  approximations  in  the  study  of  queueing 
systems  was  introduced  in  1961  by  Kingman  [15].  Subsequently  Iglehart  [11], 
Kingman  [16],  and  Newell  [24]  provided  substantial  results  in  1965.  These 
approximations  have  also  been  applied  to  problems  found  in  such  diverse 
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fields  as  statistics,  engineering,  physics,  genetics  and  neurophysiol- 
ogy. An  interesting  historical  review  of  the  use  of  the  diffusion 
equation  is  in  reference  14. 

When  employed  in  describing  queueing  systems  with  congestion  or 
heavy  traffic,  diffusion  approximations  provide  very  accurate  results. 

A system  is  considered  congested  when  the  traffic  intensity,  P , is 
never  much  less  than  unity  (say,  p > .70),  where  the  traffic  intensity 
generally  measures  the  ratio  of  the  system's  input  rate  to  its  output 
rate.  The  system  we  are  modeling  should  fit  into  this  category.  The 
fixed  cost  of  each  channel  is  extremely  high,  and  thus,  the  Imputed  cost 
of  server  idle  time  is  high.  Therefore,  it  is  reasonable  to  expect  the 
number  of  service  channels  will  be  kept  to  a minimum,  forcing  the 
traffic  intensity  to  be  high  in  many  real  situations. 

The  previous  work  done  in  approximating  the  line  length  distri- 
bution for  multi-server  queues  can  be  divided  into  three  categories: 

1)  p < 1,  2)  p > 1,  and  3)  P - 1.  The  third  category  is  highly 
unlikely  in  practice,  so  this  case  will  not  be  discussed. 

In  the  case  p < 1,  Iglehart  [11]  developed  approximations  for  the 
M/M/C  queue  and  a machine- repairman  problem.  Although,  he  provides  weak 
convergence  results  in  some  extreme  cases,  his  approximations  are 
generally  not  very  accurate  because  he  did  not  restrict  the  queue 
lengths  to  be  non-negative.  Halachmi  and  Franta  [10]  incorporated  this 
restriction  into  their  analysis  and  demonstrated  good  approximations  for 
the  GI/M/C  queue.  Their  work  will  be  discussed  in  greater  detail  later 
in  this  section.  Fischer  [6]  introduced  an  approximation  method  for  the 
distribution  of  the  virtual  waiting  time  in  an  M/M/1  queue  subject  to 
breakdowns.  His  analysis,  however,  cannot  be  extended  to  multiple 
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server  systems  nor  does  it  help  in  approximating  line  length  distri- 
butions . 

When  p > 1,  the  number  of  units  in  the  queue  will  become  unbounded, 
almost  surely,  as  time  goes  to  infinity.  Although  a stationary  prob- 
ability distribution  does  not  exist,  we  can  explore  the  transient  dis- 
tribution. The  non-stationary  problem  is  not  discussed  until  Chapter  5. 
For  continuity  of  development,  these  diffusion  approximation  methods 
will  be  covered  in  this  section.  Iglehart  and  Whitt  [12]  developed  a 
diffusion  approximation  to  the  transient  distribution  of  the  line 
length  for  a GI/G/C  queue.  Their  results  are  accurate  for  the  case  when 
p >1.  Although  they  prove  some  weak  convergence  results,  Newell  [27] 
provides  an  improvement  to  the  transient  distribution  for  all  categories 
of  traffic  intensity.  Both  of  these  approaches  will  be  covered  in  more 
detail  later. 

Diffusion  approximations  are  developed  by  essentially  replacing  the 
queueing  system's  discrete  state  space  with  a continuous  state  space. 
Conditions  are  imposed  on  the  continuous  state  space  so  that  the  newly 
defined  process  captures  the  characteristics  of  the  original  process. 

When  modeling  a queueing  process  in  this  manner,  one  gets  a partial 

/ 

differential  equation  with  boundary  conditions  that  play  a rather 
natural  role.  This  differential  equation  is  called  the  diffusion 
equation  and  its  solution  will  yield  a probability  density  function. 

Since  we  are  trying  to  find  a distribution  describing  a discrete  process, 
the  density  function  will  have  to  be  integrated  over  specific  intervals 
to  yield  the  desired  approximation. 

In  this  section  we  first  derive  the  basic  diffusion  equation  and 
the  required  boundary  conditions.  We  will  then  review  the  previous  work 
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done  in  obtaining  diffusion  approximations  to  the  line  length  distri- 
butions for  multi-server  queues.  Then  an  approximation  w ill  be  derived 
for  the  stationary  line  length  distribution  for  a multi-server  queue 
subject  to  server  breakdown.  Numerical  examples  will  be  given  to 
compare  this  approximation  to  analytic  and  simulation  results. 

3.3.1  Derivation  Of  The  Diffusion  Equation 

This  section  will  discuss  some  of  the  underlying  assumptions 
necessary  to  develop  a diffusion  approximation,  introduce  some  new 
notions  and  notation,  and  then  display  the  derivation  of  the  diffusion 
equation  and  the  boundary  conditions  imposed  on  it.  The  derivation 
follows  one  found  in  Kleinrock.  (See  [17],  pp.  69-71.) 

The  first  assumption,  which  was  previously  mentioned,  is  that  we 
will  only  be  examining  queues  with  heavy  traffic.  Thus  it  is  reasonable 
to  expect  that  the  number  of  units  in  the  system,  N(t),  is  relatively 
large  compared  to  unity.  On  a coarse  scale  of  measurement,  N(t)  changes 
very  little  in  a short  period  of  time.  Although  N(t)  is  discrete,  it  is 
mathematically  convenient  to  view  it  as  a continuous  random  variable, 
X(t),  and  thereby  allow  "infinitesimal"  queue  changes.  To  be  consistent 
with  the  nature  of  the  queue  being  modeled,  we  will  define  X(t-)  as  a 
continuous- time,  continuous-state  Markov  process  with  conditional 
transition  probability 

P(w,x;y,t)  - Pr[x(t)  _<  y | x(x)  «w]  for  T < t.  (3.46) 

So  F(w;f;y,t)  is  the  probability  that  the  process  X(t)  is  no  greater 
than  state  y at  time  t given  that  it  was  in  state  w at  time  x . 

We  will  assume  that  the  conditional  probability  density  function, 
f(w,x;y,t),  exists,  is  continuous  and  twice  differentiable.  Then  we 
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have 


_ * \ 3F(w,x;y ,t) 

f(w,T;y,t)  - * — ^ — 

3y 


(3.47) 


This  density  function  satisfies  the  Chapman-Kolmogorov  equation 

f “ 

f (w,T ;y , t)  - f (x,u;y , 

J —00 


t)  f(w,T;x,u)dx  for  T<  u < t. 

(3.48) 

Define  the  conditional  mean,  M(x,t;t),  to  be  the  expected  value  of 


the  process  X at  time  t,  given  it  was  at  x at  time  T . Define  the 
conditional  variance,  V(x,T;t),  in  the  same  manner.  Thus  we  have: 

M(x,t; t)  - E[X(t)  j X (t)  - x]  for  T < t,  (3.49) 

V(x,x;t)  - E { [X(t)  - M(x,T;t)]2  |x(T)  - x } for  T < t.  (3.50) 

Notice  that  these  moments  are  dependent  on  the  state  of  the  process  as 
well  as  the  time.  This  is  an  important  distinction  which  is  valuable  in 
approximating  multi-server  queue  distributions  (as  opposed  to  single 
server  queues).  We  shall  assume  that  these  moments  have  continuous 
derivatives  which  are  defined  as  the  infinitesimal  mean,  m(x,t),  and  the 


infinitesimal  variance,  0 (x,t).  Specifically  we  have: 


, \ 3M(x,t;f) 

m(x,t)  = kzr1 — 


0 (x,t) 


3t 

3V(x, t ;t) 


3t 


T»t 


x*t 


(3.51) 

(3.52) 


Incorporating  definition  (3.49)  and  the  fact  that  M(x,t;t)  ■ x,  we  have 


m(x,t)  * lim  t—  [M(x,t;t+At)  - M(x,t;t)] 
At->0  C 


iim  fr- 
At->0  At 


f 

J — O 


yf (x, t ;y , t+At)dy  - x 


*>  lim 

At->0 


(y-x) f (x, t ;v , t+At)dy 


(3.53) 
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Similarly, 


_2 , . x ,,  f , ^2  f (x,t ;y ,t+At)  , (3. 

a (x,t)  = lim  (y-x)  x • i 1 dv  , 

At->0  A 


54) 


and,  in  general,  we  define 


E^(x,t)  a lim 


At— >0  1 


(y-x)" 

(3.55) 


where  E (x,t)  Is  the  infinitesimal  nth  moment.  These  relations  will 
n 

become  useful  later. 

Now  we  will  derive  the  forward  diffusion  equation  by  employing  an 
expedient  analytical  technique.  An  arbitrary  integral,  I,  will  be 
defined.  Then  an  alternate  representation  using  Taylor  series  will  be 
derived.  By  taking  the  difference  of  these  two  integrals  and  by  employ- 
ing a theorem  of  integration  we  will  obtain  the  diffusion  equation. 

Consider  an  arbitrary  function  Q(y)  which  is  infinitely  differen- 
tiable and  sufficiently  bounded  so  that  the  integral 


Q(y)  dy 


(3.56) 


is  well  defined.  Using  the  definition  of  a partial  derivative  and  the 
Chapman-Kolmogorov  equation  (3.48)  we  get: 


Q(y) 


lim  ^Cw^T^^t+At)  - f(w,x;y;t_)j 


dy 


At->0 


lim  At 
At->0 


r00 

Q(y) 

—00 

_ 

f(w,t;x,t)  f (x,t;y,t+At)dx 


ay 


Q(y)  f (w,T ;y , t)  dy 


(3.57) 


Let  I ■ lim  i—  (I^-Ij) , where  1^  and  Ij  are  the  two  integrals  above. 


At->0 
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Examining  1^  alone,  we  substitute  the  Taylor  series  expansion  for  Q 
about  x and  then  interchange  order  of  integration  to  obtain 


-<»n»0 


l i.  <y-x)n  ^ 

-0  n*  dxn 


f (w;T;x,t)  f (x,t;y,t+At)dx 


dy 


00 

r „ i 

f (w,r;x,t) 

y 1 dnQ(x) 

(y-x)  f(x,f,y,t+6t)dy ] 

dx 

.CD 

_n-0  n!  dxn 

— oo 

(3.58) 

Careful  examination  of  the  integrand  reveals  that  the  inner  integral  is 
the  definition  for  the  nth  moment. 

Taking  the  limit  of  ^ we  get 

r r 


Hm  At  *1  = At  I 
At-0  At-KD  J 


Q (x)  f(w,T;x,t)dx  + 


f (w,T; 


;x,t)  t) 

n=l 


dnQ(x) 


dx 


dx 


\1^  It  l2  + l,  n! 
At-K)  n=l 


r 


f (w 


,T;x,t)  En(x,t)  dx  . 


dx 


Substituting  back  into  (3.57)  we  get 


1 « l J. 
1 ^ * 
n*l 


f(w,T;x,t)  E (x,t)  ” dx  ' 

o n dx 


(3.59) 


(3.60) 


The  nth  term  (n  - 1,2,...)  in  this  equation  can  be  integrated  by  parts 
n times  to  remove  derivatives  of  Q.  For  example,  define 

fOO 

f (w,T;x,t)  E (x,t)  d dx,  n - 1,2,...  .(3.61) 

» n dxn 
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Let 


All*) 


u - f (w,T;x, t)  E (x,t)  and  dv 

n , n 

dx 


dx. 


Then 


I » uv 
n 


vdu 


f (w,T;x,t)  En(x,t)  d "(n-lf 

dx 


d(n~1)Q(x) 


dx 


(n-1) 


-^[f(w,T;x,t)  En(x,t)] 


dx 


^ d(n~1)Q(x) 


dx 


(n-1) 


[jx  ff(w»T5x,t)  En(x,t)]J  dx  , n = 1,2,..., 


(3.62) 


since  by  previous  assumption  Q(x)  and  Its  derivatives  must  vanish  at 
+ . Thus  by  an  inductive  argument  we  have 

too  oo  n n 

I - I Q(x)  I - — [E  (x,t)  f(w,T;x,t)]  dx.  (3.63) 

I , n.  ^ n n 

1 n=l  dx 

Subtracting  this  equation  from  the  original  definition  of  I,  (3.56), 
yields 


f - I ^ fir  ‘V*-0 


dx. 
(3.64) 


By  assumption,  Q(x)  was  an  arbitrary  function.  Thus  the  second  factor 
in  the  integrand  must  be  identically  zero,  giving 

Mfr.XiXjt}  , y [E  (x , t)  f (w,T;x,t) ] . (3.65) 

dt  , n.  * n n 


It  is  reasonable  to  expect  the  conditional  density,  f,  to  be  "tightly 
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concentrated"  around  the  value  x.  Thus  we  shall  assume  that  the  third 
and  higher  infinitesimal  moments  are  negligible  (i.e.,  E^Cx.t)  - 0, 
n *3,4,...).  We  finally  get  the  second  order  partial  differential 
equation 

2 

ft  = -3x  tm(x,t)f]  + y — j [a2(x,t)f].  (3.66) 

3x4 

This  equation  is  known  as  the  forward  diffusion  equation,  the  one- 
dimensional Fokker-Plank  equation  and  the  forward  Kolmogorov  equation. 
When  the  infinitesimal  mean  and  variance  are  not  time  dependent,  the 
process  X(t),  which  this  equation  describes,  is  defined  as  a stationary 
Ornstein-Uhlenbeck  process.  When  the  infinitesimal  mean  and  variance 
are  constants,  the  process  defined  is  a Brownian  motion  or  Weiner 
process  with  drift.  Thus,  in  some  cases,  the  normal  probability  density 
or  Gaussian  function  is  a solution  to  the  diffusion  equation. 

Boundary  conditions  must  be  imposed  upon  the  density,  f,  in  order 
to  assure  a unique  and  meaningful  solution  to  the  diffusion  equation. 

The  most  natural  conditions  are  to  require  f to  be  a probability  density 
which  is  non-zero  in  the  positive  quadrant.  By  this  we  mean 

f(w,x;x,t)  » 0,  for  x < 0 and  when  t - T,  x 4 w,  (3.67) 

and 

r 

f(w,x;x,t)dx  * 1,  for  t > x . (3.68) 

J0 

As  the  process  X(t)  wanders  through  its  domain,  we  expect  it  to 
spend  very  little  time  around  its  lower  limit,  zero.  Most  of  the 
probability  mass  of  X(t)  should  in  the  tail  of  the  distribution  and  so 
it  is  natural  to  impose  a boundary  condition  on  the  diffusion  equation 
to  enhance  this  concept.  We  can  take  advantage  of  the  first  two  con- 


dltions  to  derive  another  boundary  condition,  often  termed  the  "re- 
flecting barrier"  condition.  The  first  step  is  to  integrate  the  diffu- 


sion equation. 


I*00  2 

Ifdx  * - [m(x,t)f]dx  + j [a2(x,t)f]  dx 

Jw  9x 


§£■  fdx  - -m(x,t) f + j [a2(x,t)f] 

y x-y  x=y  I 

The  nature  of  the  system  we  are  modeling  suggests  that  X(t)  will  not 
become  infinite  in  finite  time.  Thus  condition  (3.68)  and  the  con- 


(3.69) 


tinuity  of  f yield 


lim  f(w,T;X,t)  - 0, 

x-+°o 


lim  f(w,T;x,t)  =*  0,  and 

X-«o 


r a 

fdx  = 3F 
0 dt 


[1]  = 0. 


By  taking  the  limit  of  (3.69)  we  get  the  boundary  condition 
lim  j-m(x,t)f  + j [a2(x,t)f ]1  = 0. 


(3.70) 


(3.71) 


(3.72) 


(3.73) 


The  diffusion  equation  and  the  three  boundary  conditions  derived 
above,  given  in  this  general  form,  have  never  been  solved  analytically. 
If  the  infinitesimal  mean  and  variance  are  expressed  as  functions  of 
only  one  variable  or  as  constants,  then  a solution  can  be  found.  The 
remainder  of  this  section  will  display  these  solutions. 


3.3.2  Previous  Approximations  For  Queueing  Systems 

The  key  to  finding  an  accurate  diffusion  approximation  for  a 
queueing  system  is  in  finding  good  representations  for  the  infinitesimal 
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mean,  m(x,t),  and  variance,  a (x,t).  In  order  to  find  these.  It  Is 
convenient  to  think  of  the  random  variable  N(t)  as  the  difference  of  two 
random  variables,  N(r.)  * A(t)  - D(t),  where  A(t)  represents  the  arrival 
process  and  D(t)  reprtsents  the  departure  process.  The  distributions 
underlying  A(t)  and  D(t)  then  determine  whether  m(x,t)  and  a (x,t)  are 
constants  or  functions  of  one  variable.  We  will  explore  three  diffusion 
approximations  which  differ  in  infinitesimal  moments  and  traffic  in- 
tensity. 


3,3. 2.1  Stationary  Distributions  For  The  GI/M/C  Queue 

The  most  accurate  approximation  of  the  stationary  line  length 
distribution  for  the  GI/M/C  queue  was  developed  by  Halachmi  and  Franta 
[10].  For  stationary  results  they  assume  the  traffic  intensity  is  less 
than  unity  (p  < 1) . They  were  able  to  capture  the  nature  of  this  system 
by  making  the  infinitesimal  moments  depend  upon  the  state  of  the  system. 
The  methodology  used  to  arrive  at  these  moments  will  be  displayed  first 
and  then  the  solution  to  the  diffusion  equation  will  be  provided. 
Finally,  some  comparative  results  will  be  discussed. 

We  start  by  assuming  that  the  number  of  units  in  the  system  is 
continuous-valued  random  variable,  X(t) . Letting  X(t)  ■ A(t)  - D(t) , we 
can  define 

M(x,C)  - li»  E[X(t±At)  -■M1)|X(C)  -_»1 
AfH) 

Ef A(t+At)  - A(t) I X ( t ) = x] 

J1®  At 

At-K) 

E[p(t+At)  - p(t)|x(t)  = x] 

Art)  At 


(3.74) 
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The  arrival  process  is  independent  of  the  state  of  the  system  so  we  can 
drop  the  conditional  statement  in  its  expectation.  Since  we  are  in- 
terested in  obtaining  steady  state  results,  we  need  to  find  the  limit  of 
(3.74)  as  t + «,  Since  A(t)  is  a renewal  (or  counting)  process,  we  can 
use  Blackwell's  Theorem  [4]  if  we  assume  that  the  interarrival  times  are 
non- lattice  (or  non-arithmetic)  random  variables.  This  theorem  yields 

lim  E[A(t+At)  - A(t) ) =*  AAt,  (3.75) 

t-**> 


where  A is  the  interarrival  rate. 

The  departure  process  does  depend  upon  the  state  of  the  system. 
Since  the  interdeparture  times  for  each  active  server  are  independent 
and  exponentially  distributed,  we  can  take  advantage  of  the  memoryless 
property  of  this  distribution.  Let  y be  the  service  rate  of  an  active 
channel.  To  accomodate  the  continuity  assumption  for  the  system's  state 
space,  we  assume  that  the  servers  act  as  independent  infinitesimal 
units.  This  allows  the  definition 


xyA(t)  + O(At) , 0<x<c 


Pr[D(t+At)  - D (t)  > 0 | X< t)  = x]  = { 


[cyA(t)  + 0 (At) , x > C,  (3.76) 

which  together  with  Blackwell's  Theorem  gives 

f xyAt,  0 < x < C 
lim  E[D(t+At)  - D(t)|X(t)  = x]  * < 
t"*°°  [ CyAt,  x > C 

=»  min(x,C)yAt.  (3.77) 

Thus,  using  (3.74)  the  infinitesimal  mean  is 


m(x)  =A  - min(x,C)y. 


(3.78) 
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Similarly,  we  can  define  the  infinitesimal  variance 


0 (x,t)  - lim 
At-*0 


Var[X(t+At)  - X(t)jx(t)  «x] 
At 


* lim 

At-K) 


Var[A(t+At)  - A(t)] 
At 


+ lim 
At-+0 


Var [D(t+At)  - D(t) |X(t) 
At 


(3.79) 

where  we  again  assume  that  the  arrival  process  is  Independent  of  the 
state  of  the  system.  Since  we  do  not  have  the  equivalent  of  Blackwell's 
Theorem  for  the  variance  of  a renewal  process  we  need  to  employ  a trans- 
formation to  fit  the  conditions  of  a known  renewal  theorem.  Define 
AA(t)  * A(T  + A t)  - A(t),  n = [t/At],  and  AA^^  * A(i*At)  - A((i-l)*At), 
where  the  brackets,  [•],  denote  the  greatest  integer  function.  Then  by 
our  assumption  that  arrivals  are  independent,  we  have 

V - Var  f-  Y AA.,1  = Var  [AA(t)]/n.  (3.80) 

Ln  i«i 


But  assuming  A(0)  ■ 0 we  also  get 


= — 2 Var [A(t) ] (3.81) 

n 


Combining  (3.80)  and  (3.81)  yields 

Var[AA(t)]  = nV  = — Var  [ A ( t ) ] = Var  [A(t)]/[t /At ] (3.82) 

n 


Taking  limits  we  get  the  equivalent  statements 

lim  Var  [AA(t) ] = lim  < Var [A(t) ]/[ t/At ] 

» lim  Jvar[A(t)]  At/t 

t-w  | 


(3.83) 
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From  renewal  theory  (see  [28],  pp.  180)  we  have 

lim  Var  IMtll  = x3a2  § 

t-«o  C 3 


(3.84) 


where  is  the  variance  of  the  interarrival  times.  So  we  finally  get 

lim  Var  [A(t+At)  - A(t)  ] a At.  (3.85) 

t-*» 

In  a similar  manner,  we  can  find  the  variance  of  the  departure 
process,  keeping  in  mind  the  dependence  on  the  state  of  the  system. 

Since  the  service  times  for  each  active  server  are  exponentially  dis- 

2 

tributed,  we  know  the  variance  a,  of  each  time  is  the  square  of  the 

a 

2 2 

expected  service  time  (i.e.,  a“  * 1/u  ).  Therefore,  we  have 

lim  Var  [D(t+At)  - D(t) |x(t)  * x]  - min(x,C)u  At.  (3.86) 

Thus,  using  (3.79),  the  infinitesimal  variance  is 

(3.87) 

The  diffusion  equation  (3.66)  derived  in  the  previous  section 
relates  a change  in  state  to  a change  in  time.  Since  we  seek  a station- 
ary distribution,  we  must  discard  the  conditioning  on  the  initial  state 
and  take  the  limit  as  t+  ® . Thus 

3f 

lim  f (w,T;x, t)  - f (x) , and  lim  = 0 

£-KX> 

are  necessary  for  the  existence  of  a stationary  distribution.  We  then 
get  the  diffusion  equation 

.2 


a“(x)  = X^a^  min(x,C)p. 

cL 


— 2[a2(x)  f(x)]  Lm(x)  f(x)]  = 0, 

dx 


(3.88) 


and  the  associated  boundary  conditions 


ik  [°2(x)  f(x)1 


- m(0)  f (0)  =■  0, 

x»0 


(3.89) 


J 
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and 


1 

f(x)dx  = 1. 

0 


(3.90) 


The  solution  arrived  at  for  this  differential  equation  by  Halachmi  and 
Franta  [10]  is 


2 u_1 

H^[a  (x) ] exp(-2x)  , 0 < x = C 


f (x)  - < 


..  (~ 2m(C)x_| 

H exp  2 

L o(C)J 


, x > C, 


where 


u = ^ [(Aa  )2  + l], 

u a 


r 

Jo 


f(x)dx  = 1, 


and 


(3.91) 


(3.92) 

(3.93) 


H1[a2(C)  ] 


u-1 


exp(-2C) 


(3.94) 


Conditions  (3.93)  and  (3.94)  solve  precisely  for  the  constants  and 

H2 . Condition  (3.94)  requires  f to  be  continuous.  The  solution,  f,  is 

a probability  density  function.  To  get  the  approximation  for  the 

distribution  of  the  number  of  units,  N,  in  the  system,  we  define 

/•n+.  5 


r\j 

p = Pr  [N  = n] 
n 


f(x)dx  , n = 1,2,...  . (3.95) 


' n-.  5 

Certain  adjustments  must  be  made  since  the  approximation  for  p^  is  not 
well  defined.  These  are  discussed  in  another  section. 

Some  comparisons  can  be  made  between  the  diffusion  approximation 
and  the  analytical  solution  for  the  GI/M/C  queue,  but  unfortunately  only 
for  n > C.  A known  queueing  result  for  the  GI/M/C  queue  [8]  states 


r™  J™.  ^ 
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that  when  p - A/Cy  < 1,  the  stationary  probability  distribution,  pQ, 
exists  and  there  is  an  r,  0 < r < 1,  such  that 

Pn  ” Ar11  , n > C.  (3.96) 
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Thus  as  p + 1,  then  s -*■  1.  Moreover, 

[18]  points  out  that  if  we  define  p = 


(3.97)  will  become 


P„ 


-f1 


P(1  - s)s 


n-1 


for  an  M/M/1  queue,  Kobayashi 
1 - p , then  the  approximation 


n = 0 

n > 1 . (3.102) 


As  well  as  being  very  close  in  form  to  (3.100),  this  suggests  that 
convergence  is  rapid  as  p 1. 

Basically  these  convergence  results  are  due  to  the  fact  that 

2 2 2 

a = 1/X”  which  simplifies  cr“(C).  This  suggests  that  the  diffusion 
approximation  developed  here  would  also  be  accurate  for  systems  which 
have  a coefficient  of  variation  for  the  arrival  process  close  to  unity 
(i.e. , a2  /X2  - 1). 

SL 


3. 3. 2. 2 Transient  Distributions  For  The  GI/G/C  Queue 

Iglehart  and  Whitt  [12]  explored  a diffusion  approximation  for  the 
transient  line  length  distribution  for  the  GI/G/C  queue.  They  con- 
sidered the  cases  P = 1 and  P > 1.  The  latter  case,  only,  will  be 
discussed.  Their  derivations  do  not  involve  the  solution  of  the  diffu- 
sion equation;  they  are  lengthy  and  involve  much  notation,  so  only  the 
major  results  will  be  displayed.  The  analysis  involves  modifying  the 
queueing  system,  finding  an  approximation  for  the  modified  process,  and 
proving  the  same  results  are  good  for  the  unmodified  process.  Since  the 
approximation  converges  it  can  be  expressed  as  a theorem. 


Theorem  3.5.  (Iglehart  and  Whitt) 

Suppose  we  are  given  a GI/G/C  queue  with  N(0)  = 0. 
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Let  X 

U 

0 

a 

Define 

then 


represent  the  arrival  rate  of  units, 

represent  the  service  rate  of  a busy  server, 

represent  the  variance  of  interarrival  times , 

represent  the  variance  of  the  service  times. 

Y2  - X302  + C UV  • If  p - - > 1, 

a s cy 


and 


lim  Pr 

t-wo 


U(t)  - (\-p)t  „ 

_ 1 

n1'2 

1/2 

(2ir) 

rx 


exp (-y  /2)dy 


(3.103) 


This  approximation  is  a direct  solution  of  the  diffusion  equation 
without  regard  for  the  boundary  conditions  or  the  dependence  of  the  in- 
fintesimal  moments  on  the  size  of  the  queue.  If  N(t)  were  allowed  to  be 
continuous,  then  it  would  represent  a Brownian  motion  process  with  drift 
which  would  give  positive  probability  to  negative  line  lengths.  Since 
the  process  "drifts"  away  from  zero , the  approximation  becomes  more 
accurate  as  time  progresses.  If  an  initial  condition,  say  N(0)  ■ k,  is 
introduced,  where  k is  large,  then  the  approximation  would  also  be 
improved. 

Newell  [27]  points  out  that  when  some  of  the  boundary  conditions 
are  ignored,  the  solution  found  for  the  diffusion  equation  may  not  be 
unique.  For  this  case,  the  diffusion  equation  is 


3f  3f  . a2  32f 
3t  " “m  3x  2 » 2 


(3.104) 


and  the  boundary  conditions  are 

f(w,0;x,t)  ■ 0,  for  x < 0 and  t ■ 0,  x + w, 
[ f(w,0;x,t)dx  - 1,  for  t ^ 0, 

Jo 


(3.105) 

(3.106) 
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and 


-mf (w,0;0,t) 


3f (w,0,x,t) 
3x 


x-0 


= 0 . 


(3.107) 


Using  the  same  notation  defined  in  Theorem  3.5,  Newell  provides  the 
following  solution 


f (w,0;x,t) 


r 

2 

-[x-w- (A-p) t ] “ 

+2x(X-p) 

/exp 

2 

+ exp 

2 

2y  t 

L Y J 

L” 

r •>! 

r ,1 

3 

exp 

-[x-KH-(A-U)  t ] “ 1 

2 (A-y) 

exp 

-[v+v+(X-y)t]“ 

dv 

9 

2 

■> 

_ 

L 2y-t 

Y J 

X 

L 2y  t 

- 

J 

(3.108) 


This  solution  satisfies  all  the  boundary  conditions  and  is  valid  for  all 
values  of  p . Notice  Iglehart  and  Whitt's  approximation,  (3.103), 
corresponds  to  the  first  term  in  this  solution.  They  actually  provide 
a solution  to  the  diffusion  equation  (3.104)  but  do  not  meet  the  bound- 
ary conditions.  The  second  term  in  Newell's  solution  is  a correcting 
factor  so  that  the  reflecting  barrier  condition  (3.107)  is  met. 

3.3.3  Approximation  For  GI/M/C  Queue  Subject  To  Server  Breakdown 

The  method  we  use  to  find  an  approximation  for  the  line  length 
distribution  in  a Gl/M/C  queue  subject  to  server  breakdown  is  the  same 
as  the  one  discussed  in  Section  3. 3. 2.1.  A representation  is  found  for 
the  infinitesimal  moments,  the  diffusion  equation  is  solved  and  the 
derived  density  is  integrated  to  get  the  approximation  for  the  station- 
ary distribution  of  the  number  of  units  in  the  system.  In  this  section, 
we  shall  also  restrict  the  traffic  intensity,  p ■ VuC,  to  be  less  than 
one.  After  the  approximation  is  developed,  a region  in  which  to  expect 
the  best  results  is  created.  Finally,  some  numerical  examples  are  given 


u 


so  that  we  may  compare  the  approximation  to  analytic  and  simulation 


results . 


3. 3. 3.1  Derivation  Of  Infinitesimal  Moments 

We  start,  as  before,  by  assuming  that  the  number  of  units  in  the 

system  is  a continuous  random  variable,  X(t).  We  then  express  X(t)  as 

the  difference  between  the  arrival  process,  and  the  departure  process, 

that  is,  X(t)  * A(t)  - D(t).  Given  that  C(t)  is  the  number  of  operational 

service  channels,  we  can  define  the  infinitesimal  mean 

m(x  k t)  - lim  E[A(t-t'h)-D(t-Hi)-A(t)+D(t)lc(t)-k,X(t)-xl 
h->0  h 

„ lim  E[A(t+h)-A(t)  ] _ Um  E[D(t+h)-D(t)  | C(t)=»k,X(t)3x] 
h-K)  h h-K)  h 

(3.109) 

Notice  that  in  this  case,  the  departure  process  is  dependent  upon  the 
number  of  operational  service  channels  as  well  as  the  state  of  the 
system.  Using  the  same  procedure  as  described  in  Section  3. 3. 2.1, 
we  have 

lim  E[A(t  + h)  - A(t)]  -X  h, 
t-*» 

and 

lim  E[D(t  + h)  -D(t)  | C(t)  - k,  X(t)  - x]  - min(X,k)y  h (3.110) 

t-*oo 

Thus, 

ra(x,k)  - X - min(x,k)y  . (3.111) 

Notice  m(x,k)  depends  on  both  x and  k.  We  can  find  the  infinitesimal 
mean  m(x)  using  m(x,k)  as  follows: 

C 

m(x)  ■ £ m(x,k)  Pr  (k  channels  operational | x units  in  the 

k»0  system} 

C 

■ £ [X  - min(x,k)y]*  Pr  {k|x} 

k-0 


(xu)*  £ Pr{k|x}  , 0 <_  x £ 1 
k-1 


(xu)  £ Pr{k | x}  - uPr{ 1 | x}  , 1 < x < 2 

k-2 

C 2 

(xu)  £ P*{k | x)  - u £ k Pr{k|x),  2 < x <_  3 
k=3  k=l 


X - u I k Pr  {k|x} 
k=l 


, x > C 


(3.112) 


The  infinitesimal  variance,  can  be  written  as  follows: 


a2(x  k t)  = lim  Var  [A(t+h)-D(t+h)-A(t)+D(t) |c(t)-k,  X(t)-x 
h-K)  h 


lim  Var[A(t+h)  - A(t)]  + lim  Var[D(t+h)  - D(t) 1 C(t)=k,  X(t)=»x] 

h-K)  h h-K)  h 


Again  we  know  that 


lim  Var[A(t  + h)  - A(t)]  = X2  a2  h, 

t-+oo  3 


(3.113) 


(3.114) 


lim  Var[D(t  + h)  - D(t)  | C(t)  =k,  X(t)  = x]  - min(x,k)  Uh. 

t-*°°  (3.115) 


a2(x,k)  * X^a2  + ®in(x,k)  u. 


(3.116) 


Taking  expectations,  we  see  that 


2 ^2  * 
a (x)  ■ £ a (x,k)*Pr  (k  channels  operational  x units  in  the  system) 

k=0 


r 
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(3.117) 

3. 3. 3. 2 Definition  of  an  Approximation  Region 

Finding  the  conditional  distribution  Pr  {k|  x}  is  as  difficult  as 
solving  for  the  distribution  we  wish  to  approximate.  Using  Bayes' 

Theorem  to  solve  for  Pr  {k|  x},  where  the  continuous  variable  X(t)  has 
been  replaced  by  the  original  discrete  variable  N(t),  we  get 


Pr  {k|n} 


Pr  (n 1 k}  Pr  {k} 

“c — 

l Pr  {n | k}  Pr  {k} 
k-0 


Pr  {k.n} 
Pr  {nF 


(3.113) 


If  we  knew  any  of  the  probabilities  in  the  above  equation,  then  we  would 
not  need  to  use  approximations  since  solving  for  Pr  {n}  is  our  objective. 

Under  some  conditions,  however,  it  may  be  reasonable  to  use  Pr{k} 
in  place  of  Pr{k|  n}.  Looking  back  at  equations  (3.112)  and  (3.117), 
one  can  see  that  the  probabilities  of  interest  are  only  Pr{k|  n} 


i 
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when  n > C.  For  example,  in  the  two  channel  case,  we  wish  to  have 

Pr{  k [n}  = Pr  {k}  , for  k * 1,2  and  n * 3,4,... 

(3.119) 

The  stationary  probability  distribution  for  the  number  of  operation- 
al channels  is  given  by  Theorem  3.1,  equation  (3.2). 

A comparison  between  the  conditional  distribution  Pr{k|  n}  and 
Pr{k)  was  done  for  one  and  two  server  systems  using  the  results  of 
Theorems  3.2  and  3.3.  It  appears  that  the  relative  error, 

Pr{k}  - Pr{k]n},  is  a function  of  X/C  , y/n  and  X/Cy,  where 
Pr{k] n} 
c 

C * £ kPr  (k } and  the  last  term  is  the  traffic  intensity,  P . As 

k-0 

P -*•  1 the  relative  error  decreases.  As  X/£  increases  beyond  a certain 
point,  the  approximation  improves.  As  y/n  increases,  the  approximation 
gets  steadily  worse,  A proposed  region  for  approximations  which  have 
less  than  ten  percent  relative  error  is  as  follows: 

p > .75,  (3.120) 

X/C  > 1.00,  (3.121) 

y/n  < 1.00.  (3.122) 

The  comparative  data  are  displayed  in  Figures  3.3,  3.4  and  3.5. 
Figure  3.3  shows  the  results  in  a system  having  one  server  and  no 
spares.  Figure  3.4  shows  the  results  for  two  servers  and  no  spares. 
Figure  3.5  compares  results  for  the  single  server,  single  spare  system. 
Negative  errors  show  up  because  one  density  does  not  necessarily  bound 
the  other  for  all  values  of  the  parameters.  The  natural  logarithm  of 
X/£  was  used  for  the  abscissa  to  show  the  insensitivity  of  the  relative 
error  to  large  changes  in  X/C  . The  numerical  results  also  exposed 
considerable  round-off  error  problems  when  trying  to  apply  Theorems  3.2 


Pr{k 


Pr {k=2 


Figure  3.5:  Relative  Errors  for  Single  Server 
One  Spare 


1 
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An  interesting  effect  was  also  observed  when  the  relative  sizes  of 
the  server  breakdown  rate  and  repair  rate  were  increased  while  main- 
taining the  same  ratio.  Suppose  we  fix  a set  of  base  rates,  £ and  q, 
and  define  /Kq  , where  K£  and  Kq  are  the  actual  breakdown  and 

repair  rates.  Then  as  K increases,  Pr  {k}  remains  the  same  (see 
equation  (3.2));  but,  the  relative  error  decreases.  These  results  are 
displayed  for  the  single  server,  single  spare  case  in  Figure  3.6.  It 
could  be  said  that  server  breakdowns  induce  a non-stationary  service 
rate  and  that  the  "more  stationary"  the  service  process  is,  the  more 
independent  the  line  length  is  from  the  number  of  operating  channels. 

This  is  similar  to  a conjecture  by  Ross  [30]  for  a single  server  system 
with  non-stationary  arrivals.  Unfortunately,  it  is  difficult  to  quanti- 
fy  this  phenomenon  and  define  a reasonable  bound. 

RELATIVE 


Figure  3.6:  Effects  of  Relative  Size  of  K£  and  Kn 
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3. 3. 3. 3 Solution  Of  The  Diffusion  Equation 

Within  the  region  defined  by  relations  (3.120),  (3.121)  and  (3.122), 
we  shall  define 


n-1  C 

X - y £ k Pr(k)  - yx  £ Pr(k),  n-1  < x <_  n,  n-l,2...,C 
k-1  k»n 


X - y l k Pr(x} 
k-1 


, x > C 


X - xy 
n n 


, n-l<  x < n,  n-l,2,...,C 


, x>C , 


(3.123) 


where 


n-1 

X - X - y £ k Pr{k} 
n k-1 


, n-1, . . . , C+l, 


(3.124) 


yn  = u t Pr(k} 

k-n 


, n-1, . . . ,C. 


(3.125) 


Similarly, 


a2(x') 


v + xu 
n n 


where 


, n-1  < x < n,  n=l,...,C, 


, x > C, 


(3.126) 


v - XJa  + y 7 kPr{k} 
n 3 . _ 

k-1 


, n-1, . . . ,C+1, 


(3.127) 


a)  - y T Pr{k} 
n 


, n— 1,2, ... ,C. 


(3.128) 


These  moments  can  now  be  used  to  solve  the  diffusion  equation  subject  to 
the  boundary  conditions,  which  are  repeated  from  Section  3. 3. 2.1, 


~ [02(x)f(x))  [o(x)f(x)]  - 0, 

dx 


(3.88) 
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1 d_ 

2 dx 


[cr2(x)f  (x)] 


x=0 


m(0) f (0) 


0, 


(3.89) 


r 1 

f (x)dx  = 1.  (3.90) 

• 0 

2 

It  is  easy  to  establish  that  m(x)  and  a (x)  are  continuous  for 
x >0.  These  functions,  however,  are  not  differentiable  at  the  lattice 
points  {1,2,...,C}  . Thus  one  must  be  careful  to  evaluate  the  diffu- 
sion equation  only  on  the  interior  of  each  interval,  (n-1,  n) , where 
n = 1,2,...,C.  Integration  of  the  diffusion  equation  yields 


\ ^ [a2(x)f(x)]  - m(x) f (x)  = Hr  , 


(3.129) 


where  is  the  constant  of  integration  peculiar  to  the  evaluation  on 
the  interval  (n  - 1,  n),  n = 1,2,...,C,  and  Hc+^  is  for  the  interval 
(C,°°  ). 

From  the  reflecting  boundary  condition,  (3.89),  we  get  » 0.  In 
order  to  have  a proper  density  function,  we  need  lim  f(x)  * 0 and 

X-KO 

lim  d f(x)  » 0.  Thus  Hc+1  ■ 0. 
x-«o  dx 

2 

Let  mn(x),  <?n(x) , and  ^n(x)  be  the  respective  evaluations  of 

these  functions  on  the  intervals  (n  - 1,  n)  for  n - 1 C,  and  (C,00  ), 

2 

for  n ■ C + 1.  Then  by  continuity  of  m(x)  and  cr  (x)  we  have  (for 
n * 1, . . . , C) 


mn(n)  " Wn)’ 

(3.130) 

oj(n)  - <^+l  (a). 

(3.131) 

To  establish  the  composition  of  f,  we  shall  require  f to  be  continuous. 


| 
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Thus 

f (n)  * *.,(10,  n - 1,...,  C.  (3.132) 

n n-t-i 

With  the  three  continuity  relations  above  it  is  natural  to  set  « 0 
for  n ■ 2, . . . , C. 

Now  the  differential  equation  (3.129)  is  a homogeneous,  first  order 

equation  with  variable  coefficients.  Using  well  known  techniques,  the 

solution,  in  general  form,  is  found  to  be 

H (3.133) 

f (x)  » — S exp  [K  (x) ] , n-1  < x < n for  n=l,...,C  and  x > C for 

n _/ , . n — _ . . 

a (x)  n«C+l , 

n 


where  K (x) 
n 


2m(x) 

cr2(x) 


dx 


and  H is  a new 
n 


constant  of  integration.  Examining  Kfl(x)  more  closely  for  n 


1 , • • • , C , 


we  have 


KnM 


2 (X  - xy  ) 
n n 


v + XU) 
n n 


dx 


-2xy 


U) 


+ 2 


X y v 

n n n 

u)  2 

n U) 

n - 


ln(v  + xu)  ) . 
n n 


(3.134) 


For  n ■ C + 1 , 


K 


C+l 


^ 2Xc+i  . 2Xc+ix 

dx  — 

C+l 


(3.135) 


C+l 


One  can  also  show  that  K(x)  is  continuous  for  x 0.  Thus  the  final 
solution  form  for  f(x)  is 


u - 

H [v  + xu)  ] n exp[-2xy  /co  ],  n-1  < x < n,  n-l,...,C, 


n n 


n n 


f(x)  -< 


HC+1  6XP  ^2X<"J-i  x/vr4.il»  x >C* 


C+l  ' C+l 


(3.136) 
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where 


n cl) 


\ + 
n 


u v 
n n 


co 


, n=l, . . . ,C. 


(3.137) 
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It  now  remains  to  solve  for  H and  evaluate  p 

n *n 

f\, 

point  out  an  ambiguity  in  evaluating  p 


First  we  shall 
We  will  use  the  method 


'U 


n-K  5 


■>n-.5 


f(x)dx  to  compute  the  approximating  density.  This  evalu- 


ation, however,  is  not  well  suited  for  p^.  Methods  for  approximating 
will  be  discussed  in  the  next  section. 

The  constants  can  be  evaluated  by  using  the  continuity  constraint 
and  the  last  boundary  condition,  which  requires  f to  be  a proper  density 
function.  If  we  define  gR(x)  to  be  the  variable  portion  of  f^(x)  in 
equation  (3.136),  then  fn(x)  = Hngn(x)  anC*  cont^-nu^ty  requires 


Hngn(n)=Vl  8n+l  (n)  ’ n " 1 C' 


(3.138) 


We  can  solve  for  Hn  in  terms  of  H ^ using  (3.138)  recursively. 

HC  = HC+1  gC+l(C) 
gC(C) 


H I - Vn8Ctl^C)  8C(-C~1) 
. 8c(C)  gc_1(C-l) 


H = HC+lgC+l  (C)*"gn+l(n)  - H„.,R  . n «*  1 C. 

" C+l  n 


8q(C)  ...  g^ (n) 


(3.139) 


Define 


c n 


gn(x)dx 

n-.5 


and  G 


n+1 


n+.  5 


gn+i(x)dx,  n-1,. . . ,C. 
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Then 


'V/ 

P_ 


•n+.  5 

•n 

m 

f(x)dx  - 

f (x)dx  + 
n 

n-.5 

n-.5 

n+.  5 

Wx)dx 


- HF  + H G . , n - 1,...,C. 
n n n+1  n+1 


(3.140) 


Assuming  we  have  a good  approximation  for  Pq,  the  boundary  condition 
(3.90)  yields 


1 ~ Pr 


f (x)dx 
.5 


s * 

\ pn  + 

n-1 


HC+1  8C+1  ^dx 
C+.5 


I [H  F + H , - G , . ] + H_ . , 
L n n n+1  n+lJ  C+l 


n-1 


exp 

C+.5 


2X<*-i* 
L VC+1  J 


dx. 


Substituting  in  (3.139)  gives 


C C v 

1 - p.  - H.^,  J T R F + T R G + G- . . - exp 

0 c+ls  L.  n n n n C+l  2A_.,  v 

n-1  n-2  C+l 


2Wc+-5) 


V 


Finally  we  have 


C+l  J 

(3.141) 


1 - Pn 


"C+l 


R1F1  + GC+1  + ^ Rn[  W ' 2^  exp 

n—4 


C+l 


[5x^T«-5T 


C+l 


C+l  J 


(3.142) 


We  can  then  use  equation  (3.139)  to  solve  for  the  remaining  Hn,  n-l,...,C, 

f\, 

and  equation  (3.140)  to  solve  for  P , n - 1,2,.... 

n 
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3. 3. 3. 4 Approximation  For  Pg 

As  pointed  out  in  the  last  section,  we  cannot  use  diffusion  approx- 

f 5 


'v 


imations  to  estimate  p accurately.  Defining  p 


0 


f(x)dx  under- 
.5 


estimates  the  value  and  empirically  we  find  that  Pg  =*  J f(x)dx  is  also 

rn+I  ’ 5 


'V, 


not  a good  estimator.  Redefining  all  estimates  p 


f (x)dx,  n = 0,1, . 


is  also  ineffective.  This  section  will  provide  some  analytic  formulas 
which  are  good  in  specific  cases  and  then  discuss  several  alternative 
estimates.  In  all  cases,  the  estimates  have  been  empirically  found  to 
be  at  least  as  large  as  Pg.  The  cases  considered  are  briefly  described 
below. 


Case  1.  M/M/1  subject  to  breakdown,  no  spares: 


[un  - A(n  + 5)1  (n  + A + Q 


y(£  + n)  + n) 


(3.143) 


from  Theorem  3.2,  equations  (3.9)  and  (3.10). 

Case  2.  M/M/2  subject  to  breakdowns,  no  spares: 

p0  * Pgg  + P-^g  + P^p,  W^ere  C^ese  probabilities  are  defined  in  Theorem 

3.3,  equations  (3.17),  (3.18),  and  (3.19). 

Case  3.  M/M/1  subject  to  breakdowns,  one  spare: 

Pq  “ Pqq  + P-^q  + P2Q,  where  these  probabilities  are  defined  in  Theorem 

3.4,  equations  (3.43),  (3.44),  and  (3.45). 

Case  4.  M/M/C  with  no  breakdowns: 


P 


0 


C-l 

l 


n=0 


— (— ) + — (— ) (^ — 

nPy'  CPy'  xy-A 


) 


-1 


(3.144) 


This  analytic  result  was  suggested  by  Halachmi  and  Franta  [10]  as  an 
approximation  for  GI/M/C  queues  when  pg  is  very  small.  An  adjustment  to 
p,  to  reflect  slower  service  due  to  failures,  will  be  provided  later. 
Case  5.  M/G/  00  • 


PQ  - e"XE[S]  (3.145) 

where  E[S]  is  the  expected  service  time.  This  equation  is  a result  of 
Palm's  Theorem  [8]  and  is  useful  when  the  number  of  servers  is  greater 
than  two.  Here  we  assume  that  when  a channel  fails,  the  unit  being 
repaired  does  not  change  channels.  Evaluation  of  E[S]  requires  some 
analysis  which  will  be  provided  later. 

Case  6.  M/G/l: 


XE [ S ] , 


(3.146) 


where  E[S]  is  the  expected  service  time.  This  is  the  analytical  solu- 
tion for  the  M/G/l  queue,  proven  in  reference  6.  No  analytic  results 
exist  for  the  M/G/C  queue.  For  a multi-channel  queue,  we  will  assume 
that  all  units  are  served  by  a single  server  which  works  as  a "super 
server"  at  a rate  Cp  , where  C is  the  expected  number  of  operational 
channels. 


Cases  1,  2,  and  3 are  analytical  results  and  can  be  directly 

% 

applied  to  solve  for  Pg.  Cases  1,  2,  4,  5,  and  6 are  analytic  results 
for  specific  systems  which  can  be  adapted  to  provide  reasonable  approxi- 
mations. Finding  a way  to  adjust  the  service  rate  (p)  is  the  key  to 
altering  the  equations  (3.143)  through  (3.146)  to  suit  our  needs.  This 
will  be  done  by  assuming  that  all  channels  are  occupied  and  finding  the 
expected  service  rate  when  servers  are  subject  to  breakdown. 
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The  total  time  a typical  unit  occupies  a channel  can  be  expressed 


as 


S + Dx  + D2  + ...  + Dn(s), 


(3.147) 


where  S is  the  service  time,  is  the  delay  caused  by  a server  failure 
and  N(S)  is  the  number  of  channel  failures  encountered  in  a time  S. 
Given  all  channels  are  occupied  and  no  other  channels  fail  while  the 
inoperative  channel  is  being  replaced,  then  we  can  find  the  conditional 
expectation 


1 1 


i-k 


(C+L-i)n 


E[Dj  k - k]  -< 


, when  k other  channels  were  oper- 
ating just  prior  to  the  failure, 
k-0,  1, . . . , C— 1 


V 


, otherwise. 


(3.148) 


Then  using  the  results  of  Corollary  3.1  and  assuming  all  delays,  D^,  are 


independent 


C-l  , 

E[D]  “ l qfc+1  F[D | K - k] , 


k-0 


(3.149) 


where  q^  is  the  conditional  probability  that  k servers  are  operating 
just  prior  to  a failure. 

We  now  need  an  expression  for  N(S),  the  expected  number  of  failures 
in  a time  S,  which  has  an  exponential  distribution  with  rate  M . Since 
failures  occur  as  a Poisson  process  having  rate  S , 


Then 


Pr  {N(S)  - k}  - 


Me 


-Ut  e 


k! 


dt. 


E[N(S) ] -Ik  Pr{N(S)  - k} 
k-0 


(3.150) 


d 
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ye 


-ut 


dt 


k-0 


- 5 


tye~wtdt  ■ C/y. 


(3.151) 


Reference  28  provides  the  result 

E[DX  + D2  + ...  + DN(S)]-E[N(S)]  -E[D], 
So,  assuming  Independence  for  all  random  variables, 

E[T]  - E[S]  + E[N(S) ] *E[D] 


(3.152) 


S + £ E[D]  -i  [l  + 5E(d] 


Therefore  we  can  define  the  adjusted  service  rate  as 


y - 


(3.153) 


(3.154) 


1 + SE[D]  ' 

f 

We  will  apply  y in  cases  1 and  2 to  approximate  p^  for  systems 
with  one  or  two  channels  and  spare  servers.  The  modified  service  rate 
will  be  employed  in  Case  4 for  all  systems  with  over  two  channels.  For 

I 

Case  5 we  can  use  1/y  as  an  estimate  for  E[S],  Case  5 will  be  used  for 
systems  with  over  two  channels  and  Case  6 will  be  used  for  all  systems 
with  spare  servers. 

Some  empirical  results  are  compared  to  values  for  p^  estimated  by 
simulations  in  Table  1.  In  all  cases,  p^  was  overestimated.  Thus  a 
reasonable  approximation  would  be  to  select  the  smallest  value. 

3. 3. 3. 5 Comparative  Analysis 

The  diffusion  approximation  derived  in  the  preceding  sections  was 
tested  against  analytic  and  simulation  results.  The  actual  stationary 


Case 


63 


probability  distributions  were  derived  from  the  results  of  Theorems  3.2, 
3.3  and  3.4.  The  simulated  distributions  were  developed  using  a simu- 
lation routine  involving  at  least  20,000  state  changes.  The  results  are 
graphically  displayed  in  the  figures  that  follow.  Each  figure  is 
labeled  using  the  following  nomenclature: 

C ■ the  number  of  channels, 

L ■ the  number  of  spare  servers  provided, 
p - the  traffic  intensity  ■ X/cy  . 

Each  graph  includes  information  on  the  ratios,  X/£  and  y/n  . 

We  know  from  Section  3. 3. 3. 2 that  the  accuracy  of  the  approximation 
Pr  {k  |n}  = Pr  {k}  improves  as  £ and  n get  large;  therefore,  we 

should  expect  the  size  of  £ and  n to  affect  our  diffusion  approxi- 
mation. As  the  relative  size  of  the  rate  of  repair  for  servers  (n  ) 
decreases,  the  ratio  y/n  gets  large.  Examining  the  data,  we  see,  the 
approximation  suffers  as  the  ratio  y/n  increases  above  the  proposed 
bound  of  one.  On  the  other  hand,  the  approximation  is  relatively 
insensitive  to  the  ratio  X/£  . Thus,  this  diffusion  approximation  is 
most  sensitive  to  changes  in  the  server  repair  rate  (n  ) and  improves  as 
the  rate  increases  relative  to  all  other  parameters. 

Some  simulation  results  show  a slight  perturbation  due  to  auto 
correlation  effects.  Regardless,  it  appears  the  exponential  tail  of  the 
diffusion  approximation  matches  the  true  shape  of  the  distribution  in 
all  cases.  When  the  approximation  is  inaccurate,  the  tail  of  the 
approximate  distribution  is  slightly  low. 
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Figure  3.8:  C=1  , I.=0,  p=.75 


igure  3.18:  C=3,  L*0,p».75 
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Figure  3.19:  03,  L=0,  p».90 


CHAPTER  IV 


AN  OPTIMIZATION  METHOD 


This  chapter  will  derive  a method  to  find  the  optimal  allocation  of 
resources  for  the  inventory  system  we  have  modeled.  We  will  first  ex- 
amine the  backorder  function  given  the  diffusion  approximation  developed 
in  the  previous  chapter  for  the  stationary  distribution.  Since  the  dis- 
tribution is  for  the  total  number  of  units  in  the  system,  it  will  be 
modified  to  describe  a multi-item  system.  Then  using  the  exponential 
tail  of  this  distribution,  we  develop  sufficient  conditions  for  convexity 
of  the  backorder  function.  Finally,  a simple  marginal  analysis  technique 
is  given  to  find  the  solution  to  the  optimization  problem,  P. 

The  optimization  problem,  P,  we  wish  to  solve  is: 
m 

minimize  ^ Ei  I (n  - p(n|  A^y.S.n.C.L) , (4.1) 

i=l  n>s^ 

subject  to 

m 

C-Cc  + L*C  + l C±Si  < B,  (4.2) 

i»l 


and 


A < yC. 


(4.3) 


Here  C,  L,  and  s^,  i = 1,  . . . ,m,  are  the  decision  variables.  Because  of 
the  high  setup  costs  for  service  channels,  it  is  reasonable  to  assume  C 
is  small.  We  shall  show  later  that,  for  fixed  values  of  C,  P is  relative- 
ly easy  to  solve;  thus,  we  shall  enumerate  solutions  for  small  values  of 
C.  For  notational  simplicity  in  the  following  developments,  the  vari- 
able C and  parameters  y,  C.  and  ri  will  be  suppressed. 


75 


4.1  Approximating  The  Backorder  Function 


After  examining  the  objective  function,  we  can  see  that  F is  not 
separable.  In  fact,  it  is  difficult  to  explicitly  express  p(n  | X^L), 
the  distribution  of  the  number  of  units  of  type  1 in  the  system,  as  a 
function  of  L.  An  expression  for  the  expected  backorders  of  item  i, 
B^Cs^.L),  will  be  created  that  will  assist  in  the  analysis. 

Define 


B^s^  L)  - l (n  - s^  p(n  | X^L). 
n>s^ 

We  now  develop  an  expression  for  an  approximation  to  p(n  | 


(4.4) 

*±,L). 


Since  the  influx  of  each  type  unit  is  an  independent  Poisson  process, 
the  entire  input  process  is  Poisson.  The  conditional  distribution  for 
the  number  of  units  of  type  i in  the  system  is  binomial,  i.e., 

Pr  (k  units  of  type  i in  system  | n total  units  in  system) 


n < k 


n-k 


n > k. 


(4.5) 


Using  the  law  of  total  probability  gives 

k 


00  ln\lXi\ 

p(k|X1,L)  =■  I\k  \X1 

n*k ' ' 


X-X, 


n-k 


n’ 


where  pn  is  the  stationary  probability  of  n total  units 
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For  n > c,  the  diffusion  approximation,  p^,  for  pn 


previous  chapter  is 


(4.6) 

in  the  system, 
given  in  the 


’ n+.5 
f(x)  dx 
. n-.5 


5 


H 


exp(Kx)  dx 


n-.5 
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where 


^ {exp  [K(n  + .5)]  - exp  [K(n  - .5)]} 
2-  [exp  (.5K)  - exp  (— . 5K) ] exp  (Kn) 
h'  r“,  n > C, 


K - 2(X  - jjCl, 

X + yC 

a'  - ? [exp  (.5K)  - exp  (-.5K)], 


(4.7) 


(4.8) 


and 


r ■ exp  (K) . 


(4.9) 

(4.10) 

Recall  that  the  average  number  of  operating  channels,  C,  is  a function 
of  the  number  of  spares  provided,  L,  so  that  K = K(L)  and  r • r(L) . 
Unfortunately,  the  expression  for  the  few  values  of  p^  where  n < C is 
not  as  concise. 

Applying  equation  (4.7)  to  equation  (4.5)  for  k > C,  we  get 

\n-k 


p<k,vu* 


'U 

p„ 


■n\IX.  lk/x-*  \“-k 


I 

n-kW  \ 


H r 


l\  \k  00  , \ 

M i 

(X-Xt)r 

^ X J n=k\  ' 

L x -J 

n-k 


(4.11) 


Under  the  assumption  that  \ < y C,  we  have  K < 0 and  thus  r < 1.  Then 

the  factor  being  raised  to  a power  in  the  summand  is  less  than  unity  and 

we  can  apply  the  binomial  expansion  to  get 

, k 


p(k|X1,L)  = 


\ _ (X-X^r 
. “ - 


k+1 


...  .. 


ha  rv ; 

" X±r  |_X(l-r)  + X±r 


k+1 


where 


and 


a >,k_1 

Ai  bi  * 


V 


X(l-r)  + X^r  ’ 


H Xbj 

-r~r 


(4.12) 


(4.13) 


(4.14) 


Equation  (4.12)  shows  the  similarity  to  the  geometric  distribution. 
Generalizing  for  all  k 0 gives 


H r 


n-k 


% 

P„ 


where 


Aibri+di,k  ■ ki°- 


C /n\  i 1 i y 

£ yk ]\X  | \ X | Pn 
n-k 


(4.15) 


i.k 


0 < k < C 


k > C, 


(4.16) 


and  and  b^  are  defined  in  (4.13)  and  (4.14). 

We  can  now  solve  for  a general  approximation  of  B^(s^,  L)  > the 
expected  backorders  for  item  i.  Since  b^  < 1,  we  have 

OO 

B1(s1,L)  - l (n-s^  p(n|X1,L) 
n-s . 
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r 


i 


t (n-si)  p(n| X^,L)  + £ (n-si)  p(n|Xi>L) 

n»C+l 


l (n-s . ) (A  b""1  + d ) + l (n-s  ) Ab"  1 
111  i*n  n-C+1  1 1 i 

n=Si 


l (n-Si)  A^1  + I (n-Si)  di>n 


n-s 


n=s , 


A,b 


i i 


(1-b^2  n“si 


+ l (n-si)  di>n_ 


(4.17) 


We  expect  B^s.^,  L)  to  be  a monotonically  decreasing  function  of 


s^.  In  fact. 


AgBi(s,L)  = Bi(s  + 1,L)  - Bi(s,L) 


- — 2 

(l-bi)/ 


. S+1_  , s]  + £ nd  - l nd 

i ij  n=s+l  ,n  n=s  ,n 


s 


Aibi 


(l-bi) 


2 V1*  - sdi,s  < 0 


since  b.  < 1 and  d^  > 0.  It  is  not  clear  that  d.  is  decreasing  in 
i i,s  — i,s 


s for  0 < s < C;  however,  since  d * 0 for  s > C,  the  function 

IfS 


A B(s,L)  is  increasing  in  s for  all  s > C.  Therefore , -B  (s  ,L)  is 

3 11 


convex  in  s^  alone,  for  s.^  _>  C. 

The  parameter  b^,  as  a function  of  L,  has  some  properties  which 
will  be  useful  later.  These  are  easiest  to  display  after  first  ex- 


7 


80 


ploring  C,  K and  r as  functions  of  L.  The  expected  number  of  operative 
channels,  C,  is  a bounded,  strictly  increasing  function  of  L.  Further- 
more, lim  C(L)  ■ C.  For  feasibility  we  assume  there  is  an  L such  that 
L-*» 

X < yC(L).  As  a result,  K is  a bounded,  monotonically  decreasing 


function  of  L.  Specifically, 

2fX-yC(L)l  2(x-yc) 

Kro  = lim  K(L)  = lim  X + yC(L)~  X+  yC  <0, 

L-*»  L-*® 

since  X < y C.  To  show  monotonicity,  we  note  that 

AK(L)  * K(L  + 1)  - K(l) 

, X-  uC(L  + 1)]  -2 [X  -yC(L)l 

X+  yc  (L  + 1)  X+  yc(L) 

-4Xy  [cq+1)  - C(L)] 

[X+yIT(L+l)  ] [X+yU(L)  ] 

since  (i(L  + 1)  > (:(L).  Since  C(L)  increases  towards  its  upper  bound  C, 
we  expect  ^C(L)  to  be  a decreasing  function  of  L.  Then  AK(L)  would  be 
increasing  (becoming  less  negative),  and  therefore  K(L)  is  a convex  func- 


tion of  L. 


The  parameter  r * exp(K)  is  a convex  function  of  a convex  variable; 
thus,  r is  convex  in  L.  It  is  bounded  below  by  r^  ■ exp  (K^  ) which  is 


in  the  interval  (0,1).  If  we  assume,  for  b. 


a continuous  variable  we  have 


X(l-r)+Xir 


, that  L is 


dBi  » Xir  + “ Xir(-Xr  + X^r  ) 

dL  [X(l-r)  + X±r ]2 


[X(l-r)  + X r] 

i jj,  1 

since  r ■ < 0.  Thus  b^  is  strictly  decreasing  in  L.  Additionally, 


lim  b. 


x<1“r«)  + X1r0O 


lies  in  the  interval  (0,1)  and  is  the  lower  bound  for  b^. 

The  above  analysis  is  displayed  graphically  in  Figures  4.1  and  4.2. 
Figure  4.1  displays  parameters  for  a system  which  is  saturated  when  no 
spares  are  provided.  In  Figure  4.2,  the  system  is  never  saturated.  The 
variable  L is  considered  a continuous  variable  for  purposes  of  illus- 


tration 


Saturated 
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4.2  An  Algorithm  For  Determining  Unit  Stocklevels  And  Service  System 
Design 

The  steps  for  a generalized  algorithm  will  be  presented  with  an 
explanation  and  justification  of  each  given  afterwards.  The  algorithm 
is  given  in  the  context  of  solving  the  problem  graphically: 


0.  Set  C 


min 


■y 

_UJ  ’ 


where  [•]  is  the  greatest  integer  function. 


1.  Determine  lower  and  upper  bounds  for  L,  L . and  L , bv 

rr  min  max 


L . = inf  {L  > Ol-r  < C(L) } and  L = inf{L  > L , |C(L)  > .95C> 

min  — 'y  — max  — min'  — 


* * 

(s  1 9 • • • 9 s } f 

m 


2.  Solve  for  the  optimal  allocation  of  spares,  S*  - 
for  a given  budget  in  item  spares  and  a mid-range  L,  using  marginal 
analysis. 


3.  Compute  B(S*,L)  for  each  value  of  L with  S*  found  in  step  2. 

4.  Select  a new  investment  level  for  item  spares  and  go  to  step  2. 

5.  Increment  C by  one  and  go  to  step  1. 

6.  Plot  fixed  total  investment  curves  and  select  the  optimal 


solution. 

The  selection  of  Cm^n  and  Lm^n  in  the  first  two  steps  provide  the 

lowest  feasible  investment  in  service  facilities.  In  a very  congested 

system,  it  is  possible  that  L . = L since  L . would  have  to  be  large 

min  max  min 

enough  to  provide  a very  high  level  of  service  reliability. 

Sti.ps  2 and  3 fix  the  investment  in  unit  spares  and  compute  the 

expected  backorders  over  a range  of  server  spares.  Since  B^(s^,L), 

i - l,...,m,  are  convex  for  s^  > C,  then  the  total  backorder  function, 
n 

B(S,L)  - £ B^(s^,L),  is  convex  for  s.^  ^ C.  In  practice,  inventory 

inl 

systems  are  often  budgeted  so  that  all  items  are  stocked  with  at  least 
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the  mean  number  of  units  expected  in  the  system  (commonly  called  the 
"expected  pipeline  inventory").  Generally,  the  mean  of  the  distribution 
p(n]  X ,L)  is  greater  than  C;  so,  we  should  be  outside  the  region 
si  <_  C,  i “ l,...,m,  provided  the  investment  budget  is  sufficiently 
large.  Under  these  conditions,  the  optimal  values  s^*,  i * l,...,m,  can 
be  found  by  marginal  analysis  due  to  the  decreasing  incremental  returns 
to  scale.  We  shall  start  with  all  set  to  some  minimum  value  and  then 
successively  increment  the  level  of  item  i*  where 

4B1Ms1.,L)  . ^ 

Ci*  Ci 

, 4 


That  is,  we  increase  the  spare  level  of  the  item  which  provides  the 
greatest  improvement  in  the  objective  per  dollar  invested.  This  allo- 
cation can  be  followed  until  the  entire  budget  is  used  since  the 
objective  function  B(S,L)  is  strictly  decreasing  in  s^,  i = l,...,m. 
Notice  this  procedure  would  not  be  altered  if  essentialities,  E^  were 

assigned  to  each  item  and  our  objective  function  becomes 

n 

B(S,L)  - l E B (s.,L). 
i=l 

We  assume  that  for  a given  level  of  investment  in  unit  spares,  the 

optimal  values  s^*,  i - l,...,m,  do  not  change  for  different  investments 

in  server  spares.  Figures  4.1  and  especially  4.2  show  that  for  small 

changes  in  L,  b.,  does  not  change  significantly.  Moreover,  we  can 

expect  the  same  relative  changes  in  b^  for  all  i * l,...,m.  We  should, 

therefore,  not  expect  great  changes  in  the  backorder  functions  for 

various  values  of  L.  The  only  area  in  which  the  backorder  function 

changes  significantly  is  in  the  range  0 £ s^  <_  C.  This  is  most  easily 

c 

explained  by  the  presence  of  the  nuisance  parameter  £ (n  - s^)d^>n, 


AD-A064  422 


UNCLASSIFIED 


CORNELL  UNIV  ITHACA  N Y SCHOOL  OF  OPERATIONS  RESEARC— ETC  F/6  15/5 
AN  ANALYSIS  OF  RECOVERABLE  ITEM  INVENTORY  SYSTEMS  WITH  SERVICE  — EtC(U) 
NOV  78  P L KNEPELL  N00014-75-C-1172 

TR-393  NL 


86 


Note  that  we  must  include  the  cost  of  establishing  C service  channels  to 
this  to  get  the  gross  investment.  The  lowest  point  on  the  trade  off 
curve  represents  the  optimal  mix  of  investments  for  this  budget.  In 
this  illustration,  it  is  best  to  purchase  one  spare  server  (at  $100,000) 
and  allocate  $400,000  to  unit  spares  when  the  investment  that  can  be  made 
is  $500,000.  With  all  the  data  displayed  in  this  fashion,  a manager  can 
plot  several  trade-off  curves  and  observe  how  expected  backorders  vary 
with  changes  in  total  investment. 


CHAPTER  V 

NON-STAT ION ARY  ANALYSIS 

5.1  Introduction  and  Motivation 

Inventory  systems  must  often  operate  in  dynamic  environments  which 
preclude  the  use  of  stationary  planning  models.  For  example,  demand  rates 
may  cycle  much  like  traffic  during  the  day  on  a city  street;  surges  in 
demand  can  occur  when  a military  environment  goes  from  peacetime  to 
wartime;  or,  expected  demands  could  steadily  increase  as  a new  aircraft 
system  is  being  purchased.  Sometimes  inventory  problems  require  a model 
for  short  horizon  planning  only,  such  as,  equipping  an  aircraft  carrier 
for  a four  month  cruise.  In  all  cases,  the  objective  is  to  describe  the 
system's  ability  to  provide  support  through  time.  This  chapter  is 
devoted  to  modeling  the  changes  through  time  of  a recoverable  item 
inventory  system,  with  servers  subject  to  failure,  where  the  unit 
demands  may  be  non-stationary . 

In  general,  we  are  considering  a system  whose  input  and  output 
processes  are  both  non-stationary.  As  will  be  shown,  it  is  extremely 
difficult  to  describe  a multi-server  system's  transient  behavior  for  the 
case  when  both  processes  are  stationary.  We  can  remove  one  element  of 
non-stationarity  by  considering  the  service  facility  as  operating  in 
different  "service  states"  through  time.  Each  state  would  represent  a 
period  when  the  output  process  is  stationary.  The  length  of  this  period 
can  be  a random  variable. 

We  will  define  the  service  state  as  the  number  of  servers  requiring 
repair.  Thus,  the  state  space  is  {0,1,..., L,L  + 1,...,  L + C}.  For 
each  service  state,  we  want  the  transient  distribution  of  the  number  of 
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units  in  the  system.  For  example,  let  G represent  the  state  set  when 
all  service  channels  are  operating,  G - {0,1,..., L},  and  B represent  the 
state  when  one  service  channel  is  inoperative,  B ■ {L  + 1}.  The  length 
of  time  the  system  is  in  each  state  will  be  represented  by  T_  and  T„. 
Figure  5.1  shows  how  the  system  behaves  through  time. 


/ 
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Figure  5.1:  System  Performance  for  Different  Service  States 

Notice,  while  the  system  is  in  state  B,  the  service  facility  is  not  as 
efficient  and  the  number  of  units  in  the  system  drifts  up.  Let  N(t)  be 
the  number  of  units  in  the  system  at  time  t.  Define  the  time  dependent 
transition  probability, 

PB(a,b,t)  - Pr{N(0)  ■ a,  N(t)  - b |service  state  B}. 

The  time  period  we  are  concerned  with,  Tg,  is  a random  variable.  Then 
the  transition  probability  for  this  period  is 


PB(a,b)  - Pr{N(0)  - a,  N(TB)  - b} 
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PB(a,b)  - /“  PB(a,b,t)  SB(t)dt,  (5.1) 

where  S (t)  Is  the  probability  density  function  for  the  variable  T . 

D O 

Using  the  same  notation,  we  can  express  the  transition  probability  of 
going  from  b units  at  the  beginning  of  service  state  G to  c units  at  the 
end  of  this  service  state  as 

PG(b,c)  - Pr{N(Tg)  - b,  N(Tb  + Tg)  - c}. 

This  analysis  can  be  carried  further  to  describe  the  transition  from  a 
to  c through  service  states  B and  G: 

00 

PBG(a,c)  - ^PB(a,b)PG(b,c)  1 (5.2) 

It  is  apparent  that  in  order  to  model  the  system  in  this  manner,  we 
must  know  two  different  probability  functions:  (1)  the  time  dependent 
transition  distribution  for  the  number  of  units  in  the  system  given  a 
particular  service  state,  and  (2)  the  density  for  the  length  of  time 
each  service  state  exists.  The  latter  density  is  frequently  referred  to 
as  the  "passage  time"  density.  We  will  explore  each  probability  £j  ction 
in  the  context  of  the  inventory  system  we  are  modeling. 

5.2  Time  Dependent  Distribution  Analysis 

The  recoverable  item  inventory  system  will  be  modeled  as  a queueing 
system,  as  discussed  in  Chapter  3.  In  this  case,  we  are  not  concerned 
with  server  failures  since  we  need  the  line  length  distribution  when  the 
system  is  in  a particular  service  state  (i.e.,  a fixed  number  of  opera- 
tive channels).  Thus,  we  want  to  describe  the  time  dependent  behavior 
of  an  M(t)/M/k  queue,  where  M(t)  is  the  shorthand  notation  for  a non- 
stationary Poisson  input  process  and  k » 0,1,..., C.  This  section  will 


provide  the  closed  form  solutions  available  for  some  special  cases  and 
then  review  some  approximation  procedures  that  are  in  the  literature. 


5.2.1  Closed  Form  Solutions 

Consider  a system  for  the  period  when  exactly  k service  channels 
are  operational,  k - 0,1,..., C.  Without  loss  of  generality,  we  will 
assume  the  period  starts  at  t - 0.  To  describe  the  system’s  time 
dependent  behavior  during  this  period,  we  need  the  transition  proba- 
bility 

Pk(a,b,t)  - Pr{N(t)  - b | 

N(0)  * a and  exactly  k service  channels  operational} 
For  the  case  k - 0,  we  have  a counting  process  with  Poisson  arrivals. 

Thus  given  a time  dependent  arrival  rate  of  X(t),  then 
PQ(a,b,t)  - Pr{b  - a arrivals  in  time  t} 

e~m(t)[m(t)]b~a  , b > a,  (5.3) 

- (b  - a)  .' 


where 


m(t)  - /*!( t)dt. 

0 

Closed  form  expressions  for  Pk(a,b,t)  for  k ^ 1 are  very  complex. 
Saaty  [31]  provides  a general  procedure  to  find  these  distributions  for 
queues  with  stationary  Poisson  arrivals.  The  procedure  involves  establish 
ing  the  balance  equations,  obtaining  a partial  differential  equation  for 
a probability  generating  function,  and  using  integral  transforms  to 
arrive  at  a transform  representation  of  the  distribution.  To  get  the 
distribution  from  its  transform  requires  a lengthy  inversion  process. 

For  the  case  k ■ 1,  we  have 


P^.b.O  - .-(i+  U,t 


Xa-b) / z 


W2''*5”  o 


oo  iv4*2lc 

, -t  f \ V (y/2) is  the  modified  Bessel  function  of  the 

where  inW  - k»  (n+k) - 

first  kind.  Obtaining  solutions  for  k > 2 is  extremely  arduous.  The 
resulting  solutions  are  not  easy  to  apply  and  are  too  complex  to  be 
useful  for  simple  comparisons.  Thus,  we  will  consider  approximation 
methods  to  arrive  at  these  distributions. 

5.2.2  Approximation  Methods 

The  application  and  accuracy  of  approximations  proposed  in  the 
literature  depend  upon  the  stationarity  of  the  arrival  process  and  the 
traffic  intensity.  Unfortunately , there  is  a paucity  of  computational 
comparisons  among  methods.  We  will  review  the  results  for  multi-server 
queues  with  (1)  stationary  Poisson  arrivals  (M/M/k)  and  (2)  non-station- 
ary  Poisson  arrivals  (M(t)/M/k). 

Kotiah  [19]  uses  an  approximate  transform  inversion  method  to 
arrive  at  an  approximate  distribution  for  the  line  length.  He  uses  the 
integral  transform  given  by  Saaty  [31]  to  arrive  at  some  numerical 
results  for  the  case  M/M/1.  His  numerical  examples  are  accurate  for 
.5  p 4 and,  in  some  cases,  provide  bounds  for  the  actual  distri- 
bution. He  describes  a method  for  generating  a distribution  for  the 
M/M/2  queue,  but  does  not  provide  any  numerical  results. 

For  congestion  cases,  Newell's  diffusion  approximation 


(given  in  Chapter  3,  equation  (3.108)  ) is  very  accurate,  especially  as 
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p exceeds  unity.  Equation  (3.108)  can  be  applied  to  cases  where  k ^ 2; 
however,  when  p < 1,  the  accuracy  suffers  since  the  infinitesimal 
moments  are  not  functions  of  the  system  size.  This  approximation  has 
computational  advantages  over  Kotlah’s  method. 

For  the  M(t)/M/l  queue,  Moore  [21]  provides  the  most  general  re- 
sults for  a single  server  queue.  He  uses  an  interactive  procedure  on  an 
Imbedded  Markov  chain  to  find  the  line  length  distribution.  His  model 
allows  for  bulk  arrivals  and  Erlang  distributed  service  times  as  well  (so 
called,  MX(t)/EY/l  queue),  and  his  numerical  examples  are  accurate  for 
.3  <_  p(t)  <_  .9. 

A computationally  simpler  approximation  is  given  by  Pokress  [29]. 

He  compares  the  finite  server  queue  (k  > 1)  to  the  M(t)/M/°°  queue  and 
gets  very  accurate  results  when  the  range  of  p(t)  is  less  than  .8.  The 
probability  distribution  for  the  number  of  customers  in  an  infinite 
server  system,  which  is  empty  at  time  zero,  is  Poisson  with  mean 

-H(t-x), 


m(t) 


X(x)  dx. 


where  \(t)  is  the  time  dependent  arrival  rate. 

Finally,  Newell  [25,26]  proposes  a diffusion  approximation  for  the 
M(t)/M/k  queue  which  is  similar  in  form  to  equation  (3.108). 

Define  F^(x,t)  ■ Pr{number  of  units  in  the  system  at  time  t is  < x) 


Pk(0,y,t)  dy. 


Then 


v* 


.(*  ferff)  - ’‘>0- 


l 0. 


x < 0, 


where 


1 

/2U 


t 


x-m(t) 

a(t) 


dy. 


and 


ft 


m(t)  ■ A 


[ky  - X Cy) 3 dy, 

0 


a2Ct) 


B + 


f [ky  + X(y) ] dy. 


0 


(5.5) 


Newell  [26]  claims  that  for  t "large  enough",  the  constants,  A and  B, 
are  negligible.  He  suggests  that  for  k > 1,  the  accuracy  improves  when 
p(t)  < 1 and  as  t gets  large.  The  approximations  of  Newell  and  Pokress 
have  computational  advantages  over  Moore's  method. 


5.3  Passage  Time  Distributions 

The  distribution  of  transition  times  from  one  service  state  to 
another  is  developed  in  this  section.  The  flows  in  the  system  are 
displayed  in  Figure  5.2. 

I r 


Figure  5.2:  Transition  Flows  for  Number  of  Inoperative  Servers 

Since  the  state  transitions  are  only  to  adjacent  states,  we  have  a 
birth-death  process.  The  states  {0,1,..., L)  are  outlined  because  all 
channels  are  in  operating  order  when  the  system  is  in  these  service 
states.  The  work  of  Keilson,  et.  al. , [7,13,14]  assume  that  the  birth 
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and  death  processes  are  Poisson.  This  section  will  discuss  their  work 
and  provide  extensions  for  Erlang  and  deterministic  death  (server 
repair)  processes. 


5.3.1  Previous  Results 

If  we  assume  the  birth  and  death  processes  are  Poisson,  then  the 
birth  and  death  times  are  exponentially  distributed.  Define  S^(t)  as 
the  probability  density  for  the  time  in  service  state  k.  Then  we  have 

Sk(0  “ (\  + Uj^e'^V^k^  , k - 0,1,...,  L + C,  (5.6) 


where  Uq  ■ \,+c  " 0.  In  addition,  the  probability  of  a transition  from 
service  state  i to  state  j is 


/ 


xi+  vi 


j - i + 1 

j - i - 1 


\ 0 , otherwise.  (5.7) 

To  best  describe  the  service  system,  we  need  the  passage  time 

densities  for  transitions  in  the  number  of  operational  channels.  For 
the  service  state  set  fo,l, . . . ,L } we  have  C operational  channels  and  for 
service  states  k > L we  have  C + L - k operational  channels.  Keilson, 
et.  al. , describe  the  passage  time  from  a "good”  state  to  a "bad"  state. 
In  our  case,  the  GOOD  state  is  defined  as  the  service  state  fo,l,...,L) 
and  the  BAD  state  is  defined  as  the  remaining  state  set  (l  + 1,..., 

L + C }.  This  is  illustrated  in  Figure  5.3.  The  server  system  can 
wander  through  the  GOOD  state  for  some  time  before  going  to  service 
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Service  State,  N(t) 


Figure  5.3:  Transitions  from  GOOD  to  BAD  States 

state  L + 1 and  thus  enter  the  BAD  state.  For  this  reason,  the  dis- 
tribution of  the  passage  time  from  the  GOOD  to  the  BAD  state  is  not  as 
simple  in  form  as  the  distributions  previously  described. 

Define  the  time  from  the  perfect  state,  Tp,  to  be  the  passage  time 
from  service  state  0 to  L + 1 and  the  post  recoverv  failure  time,  T„, 

" —————  — g 

to  be  the  passage  time  from  first  entering  service  state  L to  first 
entering  state  L + 1.  The  time  Tp  can  be  useful  when  observing  a system 
which  starts  in  perfect  operating  condition  (e.g.,  an  aircraft  carrier 
starting  a cruise) . Let 

Sp(t)  - the  probability  density  function  for  Tp, 

S <t)  * the  probability  density  function  for  Tg, 

Fp(t)  ■ the  survival  function  of  Tp 

(O0 

- Sp(y)dy, 

• t 

F (t)  ■ the  survival  function  of  T_ 

^ _ Kj 

<00 

- Sg(y)dy,  and 
t 
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V 

S.  (z)  ■ the  probe'  ‘ lity  density  function  of  the  first  passage  time  from 


service  state  k to  k + 1. 
Then  we  have 


Sp(t)  - Sq  * S+  S*(t), 


(5.8) 


and 


SG(t)  - s£(t). 


(5.9) 


where  "*"  denotes  convolution.  Using  equations  (5.6)  and  (5.7),  we  have 

(5.10) 


sj(t)  -XQ  e-X0C, 


and 


S*(t)  " xke~(Ak+yk)t  + Uke"(Xk+lJk)  *sk-l(t)*Sk(t)’  k > 1. 

(5.11) 

Graves  and  Keilson  [7]  take  advantage  of  the  convolution  properties 
of  Laplace  transforms  to  arrive  at  the  desired  densities.  Define 


O^(s)  * Laplace  transform  of  s£(t) 
Then  we  get  the  relations 


” <fsc  sj(t>  dt. 

■ 0 


x0 


0 s+x. 


<(S) 


*k 


s + \ + \ - Vk-i(s) 


Ws) 

pk+l<s) 


(5.12) 


and 


op(s)  = ir  oZ  (s), 

F k-0  K 


(5.13) 
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where  Pq(s)  ■ 1 

P^s)  - s + XQ 

Pk+l(S)  " (S+Xk+Uk)  Pk(S)  " ykXk-lPk-l(S)  * k“2,3"‘  • 

(5.14; 

Examining  the  polynomials  Pk(s)  we  get  the  following  properties: 

Theorem  5.1.  (Graves  and  Keilson) 

(i)  Pk(s)  has  k simple  roots,  qlf...fq,, 

(ii)  pk+1Cs)  has  k + 1 simple  roots,  ri» • • • >rk+1»  such  that 


< rk+l< 


qk  ' rk  " 


*1  < rl  < 


0. 


Since  pk(s)  is  a polynomial  of  degree  k,  Theorem  5.1  provides  the  com- 
plete factorization  of  P^(s).  Property  (ii)  of  the  theorem  proves  that 
the  roots  of  each  polynomial  provide  upper  and  lower  bounds  for  all  but 
one  of  the  roots  of  the  succeeding  polynomial.  Thus,  the  roots  for  each 
polynomial  can  be  determined  using  a recursive  algorithm.  From  equation 
(5.13)  we  have 


vs> 


L + 

TT  <\  (S) 
k=0 


ir  WS) 

k-0  Pk+l(S) 


XQX1’ ' ,XL 
PL+1(S) 


X0X1,,'XL 

L+l 

T T (S-r.  ) 
k=l 


(5.15) 


where  r^,  i ■ 1,.,.,  L + 1,  are  the  distinct  negative  roots  of  PL+^(s). 

Using  partial  fractions  and  the  fact  that  s~~- 
rkt  " k 

e , we  get  the  following  result: 


is  the  transform  of 
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Theorem  5.2. 


The  probability  density  function  of  the  passage  time  from  the 


perfect  state,  T , is 

P L+l  r,  t 

S (t)  - l 8,  e * , t > 0, 

P k-1  * 


where  r,  < 0 and  0.  - T ,, 
k k L+l 


X0X1*,,XL 


TT  (r.-r  ) 

i-1 

i?*k 


(5.16) 


, k-1,..., L+l 


L+l 


Notice  T 0 ■ 1;  however  8 < 0 is  possible  for  some  k.  In  a similar 


k-1 

fashion  we  get: 


Theorem  5.3. 


V 18 


The  probability  density  function  of  the  post  recovery  failure  time, 

(5.17) 


L+i  rkt 
e 


° k-1 


where  q.  , k - 1,...,L,  are  the  distinct  negative  roots  of  P (s) , 


ru,  k-1,...,  L+l  are  the  distinct  negative  roots  of  pL+1(s)» 


and 


\ ■ nr 


Wv,  W 

i-1 


L+l 


TT  (r.-r  ) 

i-1 

il*k 


In  this  case,  £ ■ 1 and  0 due  to  property  (ii)  of  Theorem  5.1. 


Notice  both  distributions  are  mixed  exponential  distributions.  This 
leads  to  some  interesting  properties  of  these  distributions.  Below  are 
some  useful  properties  given  by  Kellson  [13]. 
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Def inition 

A function  f is  log  concave  (convex)  if  ln(f)  is  a concave  (convex) 
function. 


Theorem  5.4. 

If  S^(t)  and  S0(t)  are.  log  convex  probability  densities  defined  on 

2 

some  connected  interval  T,  with  means  and  variances  o^,  then 


(i) 

Wt) 

+ X„S9(t)  is  log  convex  on  1 

> 0, 

r , x 

(ii) 

Fx(t) 

= j c S^(y)dy  is  log  convex, 

(iii) 

< 1 , and 

0* 

there 

(iv) 

is  a distribution  function  G 

S, (t)  = ye  dG(y) 
J0 


Theorem  5.5 


If  S1 (t)  and  S9(t)  are  log  concave  probability  densities  defined  on 
1 2 
some  connected  interval  T,  with  means  and  variances  a^,  then 

(i)  S^(t)  * S2(t)  is  log  concave  on  T, 

(il)  F^(t)  = S^(y)dy  is  log  concave, 

(iii)  — - 1 1 ,C  and 

afSl(t) 

(iv)  lira  — - = 0 for  some  u >0. 

t-*»  e ^ 

These  theorems  have  obvious  extensions  to  include  discrete  distribu- 
tions. The  exponential  distribution  is  the  only  function  which  is  log 
concave  and  log  convex.  Most  probability  densities  which  are  unimodal 
are  log  concave,  such  as:  the  normal,  binomial,  Poisrsn,  geometric, 
negative  binomial,  beta  and  Erlang  distributions.  The  gamma  density 
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r -x 

with  -1  < r < 0,  however,  is  log  convex. 

Theorem  5.4  proves  Chat  S^Ct)  and  TG(t)  are  log  convex  and  Theorem 
5.5  proves  that  S^(t)  and  Fp(t)  are  log  concave.  Moreover,  from  prop- 
erty (iv)  of  both  theorems,  we  can  expect  the  tails  of  S_(t)  and  S_(t) 

r 

to  be  bounded  by  negative  exponential  curves.  Since  both  S_ Ct)  and 

yj 

Sp(t)  are  mixed  exponentials,  we  expect  F"G(t)  and  F^(t)  to  have  expon- 
ential tails.  The  numerical  examples  of  Graves  and  Keilson  [7]  verify 
this.  A pair  of  typical  survival  functions  are  illustrated  in  Figure 
5.4. 


Figure  5.4:  Log  Survival  Function  Comparison  [7] 

5.3.2  Extension  to  Erlang  Distributed  Repair  Times 

We  now  consider  a system  with  a Poisson  failure  process  and  Erlang 
distributed  repair  times.  A failed  server  will  require  k "phases"  of 
repair,  each  exponentially  distributed  in  length.  For  this  case,  the 
service  state  space  must  be  modified  so  that  the  arrival  and  departure 
processes  can  be  accurately  expressed  as  functions  of  the  state  of  the 
system.  We  will  use  the  method  of  phases  and  will  assume  that  the 
server  repair  facility  has  only  one  repairman  whose  repair  rate  is 
proportional  to  the  number  of  units  in  the  system.  It  is  sufficient  to 
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expand  the  service  state  space  to  the  number  of  phases  in  the  system, 
{0,1, . . . ,k,k+l, . . . ,2k, • . . ,k(L+C) }.  Now  the  GOOD  state  is  the  set  G - 
{0, 1, . . . ,kL }and  the  BAD  state  is  the  set  3 = {kL+1, . . . ,k(L+C) }.  The  3AD 
state  can  only  be  entered  upon  the  arrival  of  a failed  server  so  we  will 
examine  the  state  of  the  system  at  each  arrival.  After  describing  the 
state  changes  between  arrivals,  we  will  construct  a transition  matrix. 
Using  the  structure  of  this  matrix,  we  can  derive  the  passage  time 
distribution. 

The  service  state  changes  at  each  arrival  epoch,  X^  , are  illus- 
trated in  Figure  5.5.  Let  N(t)  be  the  service  state  at  time  t.  Notice 


that  N(X^)  ilk,  for  all  X^  , since  the  state  is  examined  just  after  the 
arrival  of  k phases.  We  know  the  interarrival  times,  X^  - X^.  are 
independent  and  exponentially  distributed,  with  rate  C5,  while  the 
system  is  in  the  GOOD  state.  The  number  of  phase  completions  in  an 
interarrival  period  is  dependent  upon  the  initial  service  state.  The 


k(L+C) 


1 2 
_ T 


xxx 

3 4 5 


Figure  5.5:  Service  State  Changes  Between  Arrival  Epochs 
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Incerdeparture  time  for  a system  in  service  state  i is  exponentially 
distributed  with  rate  [i/k] *n,  where  [i/k]  represents  the  number  of 
servers  being  repaired.  Thus,  the  next  service  state  transition  is 
governed  by  the  following  probabilities: 

d^  ■ Pr  {a  phase  completion  is  the  next  event  ! in  service  state  n} 


[n/k]n 

C£+[n/k]n 


, 0 < n < kL 


[n/k]n 

(C+L-[n/k]£+[n/k]n 


, kL+1  < n < k(L+C) , 


(5.13) 


U ■ Pr  { a server  failure  is  the  next  event  j in  service  state  n) 
n 


- 1 - d . 
n 


(5.19) 


Each  service  state  transition  is  independent  of  the  previous  one  because 
of  the  memoryless  property  of  the  exponential  distribution.  At  an 
arrival  epoch,  we  add  k phases  to  the  current  service  state.  Thus  we 
have  the  following  transition  probabilities: 

P^j  ■ Pr  {in  service  state  j at  an  arrival  epoch  j in  service  state  i 
at  the  previous  epoch} 

“Pr  {1+k-j  phase  completions  in  an  interarrival  period, 
i + k - J ^ 0 ] start  in  service  state  i) 

- d±  •*J_fcfl«J_k 


i 

'TT 

n-j-k+1 


d U.  . , 
n j-k 


i > j - k. 


(5.20) 
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The  transition  probability  matrix  is: 


0 * • ' 0 P0k  0 


0 ...  0 


0 ...  0 Plk  pljk+1  o . . . 0 


0 • ' • 0 PkL,k PkL,kL 


0 ...  0 P 


kL+1 ,k  . . . . P 


kL+l,kL 


0 • • • 0 Pk(L+C),k  ' ‘ 'Pk(L+C),kL 


PkL,kL+l  - ' * PkL,k(L+C) 


P kL+1 , kL+1  • ' PkL+l,k(L+C) 


Pk(L+C) ,kL+l"  ,Pk(L+C) ,k(L+C; 


ch  n 

a^  represent  the  l row  of  An;  n = 0,1,..., 

b^  represent  the  i row  of  B,  i = 0,1,..., kL,  and 

1 represent  the  kC  - dimension  column  vector  (1,1,. 


Notice  that  submatrix  A represents  transitions  from  the  GOOD  state  to 
the  GOOD  state  and  submatrix  B represents  transitions  from  the  GOOD 
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state  to  the  BAD  state.  If  we  define  N to  be  the  number  of  arrivals 

P 

required  to  enter  the  BAD  state  from  the  perfect  state  (e.g.  = 4 

in  Figure  5.5),  then 


Pr{Np  - 1} 


Pr{Np  - 2} 


k(L+C) 

7 P - V p -h*l 
j£B  UJ  j-kL+1  3 U 

kL  k(L+C) 

l l P0i  Pij  ‘l  P0i  l Pij 

i£G  j£B  U 1J  1=0  U1j=kL+l  J 

kL 

j0  P01<bl'i)  ■ a0  B — ’ 


(5.21) 


Pr(N  - 3}  = l l l P P P 
ieG  j£G  k£B  J J 


P0i<L  Pi1  V “ l P0i(aL’B‘i) 


ieG  j£G  keB  i£G 


" aQ  B‘I* 


and,  in  general, 


Pr{Np  - n)  = aQn_1) #B*1, 


n = 2,3,.... 


(5.22) 


Similarly,  defining  as  the  number  of  arrivals  required  to  enter  the 
BAD  state  just  after  entering  the  GOOD  state,  we  get: 


Pr{N 


G “ 1}  " j£B  Pkl*J 


bkL  ’ 


(5.23) 


Pr{NG  = 2} 


I I pkLiPij 
i£G  j£B  ’ 

k (L+C)  kL 

l P.T  . ) - l PV1  , (b  *1) 

ieG  kL,ij=kL+l  i;1  i=0  kl,i  1 


akL  B*-’ 


■ 3t  ;L  jo  iBp^ ^ 


l pkl  id  pi  1 l pi  k> 

ieG  j£G  1* ^keB 


and,  in  general 


IV.  i-. 

I pu 

i«0  x 


akL  'B'i* 


Pr(Nc  = n>  =*  a^  L)  *B*1,  n = 2,3, 


(5.24) 


The  distribution  for  T and  T„  can  be  derived  from  the  above  work. 

P G 

For  example,  if  Np  = n, then  the  time  from  the  perfect  state,  Tp,  is  the 
sum  of  n identically  distributed  exponential  random  variables,  i.e.,  Tp 
is  Erlang  distributed  with  mean  n/(£  and  shape  parameter  n.  Thus,  we 


F (c)  - Pr(T  > t> 
P P 


| 
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£ 

- 1 - l Y (t)  Pr{N  - n},  t > 0, 

n-1  n p 


where 


V‘>  - 


2 n 


tn"1exp(-n2t/COdt. 


Similarly, 


F (t)  - 1 - l Y (t)  Pt{N  - n},  t > 0 


(5.25) 


(5.26) 


(5.27) 


n=l 

where  Y (t)  is  defined  above, 
n 


5.3.3  Extension  To  Deterministic  Repair  Times 

If  we  assume  the  repair  times  are  fixed  in  length,  the  queueing 
process  becomes  non-Markovian;  thus,  describing  the  service  state  of  the 
system  through  time  is  an  arduous  task.  The  problem  can  be  simplified 
by  considering  the  state  of  the  system  at  certain  lattice  points  in 
time  - specifically,  at  integer  multiples  of  the  repair  time,  R.  In 
this  case  the  GOOD  state  is  the  set  G - {0,1,..., L}  and  the  BAD  state  is 
the  set  B ■ (L  + 1,...,  L + C}.  The  state  transitions  and  passage  times 
are  illustrated  in  Figure  5.6.  We  will  assume  that  repair  is  only 
initiated  on  the  failed  servers  present  at  the  beginning  of  the  repair 
period.  This  could  be  viewed  as  a periodic  review  inventory  model  where 
orders  are  placed  for  new  servers  every  R days  and  the  lead  time  is  R 
days.  Notice  the  BAD  state  can  be  entered  during  a repair  period  and 
exited  at  the  completion  of  the  period.  As  a result,  we  will  have  to 
examine  the  state  of  the  system  just  prior  to  the  end  of  a repair  period. 

We  will  use  a method  similar  to  the  one  used  for  Erlang  distributed 
repairs.  A transition  probability  matrix,  P,  will  be  created,  in  order 
to  find  a distribution  for  the  number  of  service  periods  in  a passage 


time.  Then  the  results  will  be  refined  to  allow  for  transitions  into 
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Figure  5.6:  Service  State  Changes  3etween  Service  Epochs 


the  BAD  state  at  non-lattice  tines.  Define 

?ij  serv^ce  state  j just  prior  to  completion  of 

repair  period  [ in  service  state  i at  beginning 
of  repair  period  } 

* Pr  {(j-i)  servers  fail  in  time  R| in  service  state  i 
at  beginning  of  service  period  :. 

As  long  as  the  system  remains  in  the  GOOD  state,  the  servers  will  fail 
as  Poisson  events  with  rate  C Thus,  ,e  number  of  failures  in  time  R, 
given  j < L,  is  Poisson  distributed,  i.e., 


e~CCR(C£R) 

(j-i)! 


j-i 


0 < i < j < L. 


(5. 23) 

Since  the  server  failure  rate  is  constant  when  j above  is  in  Che  GOOD 
state,  equation  (5.28)  is  the  probability  of  (j-i)  server  failures  in 
time  R.  Thus,  we  have  the  following  important  property: 

P^  = Pq  Qj_i)  whenever  j ■ 0,1,..., L and  i£  j.  (5.29) 
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The  remaining  probabilities,  however,  are  not  as  simple  because  while  the 
system  is  in  the  BAD  state  the  server  failure  rates  are  dependent  on  the 
service  state. 

The  interarrival  times  are  independent;  therefore  the  number  of 
server  failures,  N(t),  is  a renewal  process.  Define  the  distributions 
Fin(t)  » Pr(the  nth  failure  occurs  at  time  £ t| 

i servers  are  being  repaired  at  time  t =*  0), 


and 


Fj(t) 


Pr(the  interarrival  time  between  failure  j and 
failure  (j+1)  is  £ t),  j - 0,1..., C+L. 


Then 


F+  * F+  * 

1 i+1 


••  * Vi^’ 


(5.30) 


where  * denotes  convolution. 

For  the  one-step  transition  distributions,  we  know 


Fj(t) 


1 - e-C , 0 < j < L 


1 _ e"(C+L  L+1  < j < L+c 


(5.31) 


Using  the  Laplace  transforms  of  equations  (5.30)  and  (5.31)  we  can  find 
the  distributions  Fij(c)i  J “ + 1,...,  L + C. 

From  renewal  theory  [28]  we  now  have  the  remaining  probabilities 
Pj,j  = Pr{(j-i)  server  failures  in  time  R| 

i servers  are  being  repaired  at  the  beginning 
of  the  service  period} 


lpii(E)  - * 


L,  L+1, . . . , C+L-l, 


j = L + C. 


(5.32) 
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Now  a transition  probability  matrix,  P,  can  be  established: 


D 

koo 

P01  • 

• * P0,L 

1 

1 

I 

P 

0.L+1 

. p 

0 ,L+C 

0 

• 

P11 

P1,L 

1 

1 

1 

1 

p 

1,L+1  . . 

. P 1 

l.L+C  | 

j 

I . 

0 

• 

‘ ° PL,L 

I 

1 

1 

1 ' 

p 

L,L+1 

p 

L , L+C 

0 

• 

. . 0 

1 

1 

1 

PL+L,L+1  * 

‘ PL+1,L+C 
• 

0 

• 

. . 0 

1 

1 

1 

1 

0 . . 

1 

* p L+C , L+C  ; 

1 

1 

1 

_! 

Let  represent  the  i row  of  the  submatrix  A, 

b^  represent  the  iC^  row  of  the  subtnatrix  B,  and 

T 

^represent  the  C-dimensional  column  vector  (1,1,..., 1)  . 

Let  Np  be  the  number  of  repair  periods  prior  to  entering  the  BAD  state 
given  the  initial  service  state  was  zero  (the  "perfect"  state).  For 
example  Np  = k in  Figure  5.6. 
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Define  Che  matrix 


p 

. . p„ 

p 

p 

00 

01 

0,L-2 

0,L-1 

01 

00 

P01  • 

‘ ' P0,L-2 

p 

0,L-1 

0 

00 

P01  • 

p 

- 0,L-2 

0 

0 

# 

• 

• 

00 

P01 

• 

0 

00 

0 . 

. 

0 

0 

where  w^  is  Che  row  of  W.  Then  using  the  property  (5,29),  we  have 


pr(Sp  . o>  - JB  p01 


L+C 

J+1  P»‘ ' V 1 


(5.33) 


L+C 


L 

ieG 

L 

j£B 

01  ij 

4 1 Oi 

1=0  U1  j=l 

L 

y 

P„, 

b . *1  = a 

*B*1  . 

L 

1=0 

Oi 

l — 

0 - 

I 

ieG 

I 

jeG 

l p0i 

keB 

P.  . P.  . . 
ij  J-i,k 

L 

L 

L+C 

y 

p„ 

y p 

y p 

L 

1=0 

01 

4 ri-j 

J-i  J 

k=L+l  J_i'k 

L 

L-i 

L+C 

I 

1=0 

P0i 

l P0  i 

j=0 

l Pjk 

k=L+l  2 

L 

L-i 

L 

y 

P„, 

y prtj 

b *i  = y p 

L 

1=0 

01 

jio  °J 

1 — > r.  1 

J i=0 

ao 

•w*B 

•1 

and  by  an  inductive  argument,  we  have 
Pr{N  ■ k)  ■ an*wk  ^’B‘1, 

p 0 — 


01  i 


k = 1,2, .. . 


(5.34) 


Ill 


Let  N^  be  the  number  of  repair  periods  from  the  time  the  system 
"recovers"  from  a BAD  service  state,  to  just  prior  to  re-entry  into  the 
BAD  state.  For  example,  in  Figure  5.6,  * n-m.  The  service  state  at 

the  beginning  of  our  observation  is  unknown;  it  is  the  number  of  server 
failures  during  the  last  repair  period  which  contained  a BAD  service 
state.  This  situation  is  expanded  in  Figure  5.7. 


Figure  5.7:  GOOD  to  BAD  to  GOOD  Stace  Transition 


Since  the  server  failure  rate  is  dependent,  the  state  j depends  upon  the 
previous  state  i.  (Notice,  service  state  i need  not  be  in  G.)  We  are 
assuming  that  the  number  of  servers  in  repair  at  time  (m-l)R  does  not 
depend  on  any  initial  conditions  (i.e.,  we  assume  steady  state  condi- 
tion). ThenfPr  N[(m-1)R]  = i}  *TT^,  where  fp  is  the  stationary  proba- 
bility of  being  in  service  state  i.  We  want  the  probability 
■ Pr  (N(mR)  = j,  jeG,  N[(m-1)R]  = i,  i + jsB} 


'j 


* £ Pr{j  server  failures  in  time  R|  N[(m-1)R]  = i}.TT. 

i+j  eB  1 

L+C-j 


i-L-j 


(5.35) 
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Let  v - (v, , . . . ,vT ) 1 . Then  in  a manner  similar  to  that  used  for  N , we 
— J.  L p 

v * B*l, 


have  L L+C 

Pr{N  = 0}  - l l v P 

i-0  j-L+1 


ij  - 

Pr{NG  - 1}  - l l l vt  P,.  P. 


i£G  jeG  k£B 


ij  (J-i),k 


L L L+C 

i-0  j-i  J k=L+l  U 


L-i 


L+C 


l vi  l P0j  l P1  k 

i-0  j-0  J k-L+1 


v*w*B*l, 


and,  in  general, 

Pr(NG  - k>  - v * wk*B’l,  k - 0,1 (5.36) 

We  now  have  distributions  on  the  number  of  whole  repair  periods  in 

the  times  T and  T_.  It  remains  to  find  the  distribution  for  the  time 
P G 

within  a period  until  the  BAD  state  is  entered.  This  time  is  repre- 
sented in  Figure  5.7  as  T„_.  This  passage  time  is  dependent  upon  the 
service  state,  i,  at  the  beginning  of  the  period,  where  i£  G.  We  have 
the  stationary  probability  „ Pr{N[(m_1)R]  „ ±\  ieG} 

_ , i-0,1, . . . ,L.  (5.37) 

L 

l i 

i-0 

The  times  in  each  state  i,  i + 1,...,  L are  independent  and  exponen- 
tially distributed  with  mean  1/C£.  Thus  the  passage  time  distribution 
is 

Pr  { TGB  < t | N[  (m-l)R]  - i } - \_i+1(t),  i - 0,1,..., L, 


(5.38) 
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where  Y^(t)  is  the  Erlang  distribution  with  mean  k/CF  and  shape  para- 
meter k.  Combining  (5.37)  and  (5.38)  yields 


Pr  ( Tgb  < c>  = l lriYL-i+l(t)' 

GB  - i=0  1 L 1+1  (5.39) 

Combining  equations  (5.34),  (5.36)  and  (5.39)  gives  us  the  desired 
distributions  (for  k = 0,1,...) 

Pr  { T < t } = Pr  { N = k}*Pr  < t},  kR  < t < (k+l)R, 


P - 


and 


(5.40) 


Pr  ( Tg  £ t}  = Pr  (Ng=  k}*  Pr  {Tgb  < tl  kR  £ t £ (k+l)R. 


(5.41) 

It  may  be  more  realistic  to  allow  repair  to  be  initiated  on  a 
server  as  soon  as  it  fails  as  opposed  to  waiting  until  the  end  of  the 
repair  period.  Then  conceptually,  if  we  consider  the  system  at  lattice 
points  only  (every  R time  units),  the  analysis  could  be  very  simple. 
Suppose  PG  is  the  probability  of  being  in  the  GOOD  state  at  the  end  of  a 
repair  period,  independent  of  the  starting  service  state.  Then  the 
distribution  for  N„_ , the  number  of  repair  periods  in  the  GOOD  state 
before  reentering  the  BAD  state,  would  be 


Pr  { Ngb  - k}  = PG  (1-PG),  k = 0,1,....  (5.42) 

Unfortunately,  PG  is  dependent  upon  the  starting  state  which  complicates 
the  analysis  considerably. 

In  the  previous  analysis,  the  service  state  transitions  could  only 
be  in  one  direction  between  lattice  points.  If  repair  on  a failed 
server  begins  immediately  and  not  at  lattice  points,  then  the  service 
state  can  wander  in  either  direction.  Figure  5.8  shows  how  the  system 
can  conceivably  enter  and  exit  the  BAD  state  during  a period  of  length 
R.  This  transition  would  change  the  server  failure  rate  and  thus  affect 


Figure  5.8:  Transitions  When  Repairs  Initiated  Immediately 

the  number  of  server  failures  in  the  period  R.  Since  the  repair  time  is 
R,  the  number  of  servers  in  repair  at  the  end  of  the  period  is  precisely 
the  number  that  failed  during  the  period.  To  find  the  distribution  for 
this  number,  we  need  to  know  (1)  how  many  servers  are  being  repaired  at 
the  beginning  of  the  period  and  (2)  how  much  repair  remains  to  be  done 
on  each  server.  Given  certain  considerations,  these  problems  may  not  be 
significant  factors. 

Suppose  that  the  service  channel  failure  rate  is  small  compared  to 

the  repair  rate  and  the  chance  of  entering  and  exiting  the  BAD  state  in 

a single  period  is  negligible.  Define  P^,  as  the  probability  of  staying 

in  the  GOOD  state  in  a time  interval  R.  If  N(t)  is  the  number  of  server 

failures  in  time  t while  in  the  GOOD  state,  then 

L 

P - Pr  (N(R)  < L}  - l Pr  {N(R)  »i}  (5.43) 

G i-0 

Since  we  remain  in  the  GOOD  state  for  the  entire  period,  then  N(R)  is 
Poisson  distributed  (equation  (5.28))  and  we  can  apply  equation  (5.42) 
to  get  a simple  result. 

We  have  now  completely  described  the  passage  time  distributions  for 
the  service  states  given  exponential.  Erlang,  and  deterministic  distrib- 
uted repair  times.  These  distributions  can  be  integrated  with  actual  or 


approximate  time  dependent  line  length  distributions  to  describe  the 
inventory  system's  time  dependent  behavior. 


CHAPTER  VI 


CONCLUDING  REMARKS 


This  thesis  presented  stationary  and  non-stationary  (time  depend- 
ent) analyses  for  a recoverable  item  inventory  system.  It  differs  from 
the  previous  work  done  in  this  area  because  the  capacity  of  the  service 
facility  was  limited  and  variable  due  to  server  breakdowns. 

The  mathematical  modeling  of  the  system  in  Chapter  II  demonstrated 
the  need  to  derive  a representation  for  the  distribution  of  the  number 
of  units  in  the  system.  In  Chapter  III  we  showed  that  analytic  attempts 
to  find  this  distribution  yield  neither  comprehensible  closed  form 
solutions  nor  a means  for  comparison  between  systems.  The  development 
of  a diffusion  approximation,  however,  did  give  us  some  insight  into  the 
nature  of  the  system's  performance.  The  approximate  distribution  was 
found  to  have  a geometric  tail  which  led  to  a simple  way  to  obtain  a 
solution  to  an  optimization  problem.  The  solution  of  this  problem  can 
be  used  by  managers  to  allocate  limited  resources  optimally  among  item 
spares  and  repair  facilities. 

The  time  dependent  analysis  of  Chapter  V provides  a basis  for 
examining  the  short  term  or  non-stationary  behavior  of  the  recoverable 
item  inventory  system.  A study  of  the  distribution  for  the  times  be- 
tween service  channel  failures  gives  a technique  to  analyze  the  effects 
of  different  service  facility  designs.  For  example,  the  frequency  and 
duration  of  changes  in  the  service  facility's  performance  can  be  cal- 
culated. This  information  can  also  be  integrated  with  transient  line 
length  distributions  to  provide  a measure  of  overall  system  performance 
through  time. 
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There  are  several  different  areas  where  this  work  can  be  extended. 
One  obvious  extension  would  be  to  consider  multi-echelon  systems.  The 
result  would  be  a modification  of  Sherbrooke's  METRIC  model  [32]  to 
include  a finite  number  of  servers.  In  the  same  manner,  items  which 
require  service  could  be  composed  of  sub-units  which  need  to  be  re- 
placed, repaired,  and  stocked.  This  multi- indentured  concept  is  con- 
sidered for  the  adequate  server  case  in  Muckstadt's  MOD-METRIC  model 
[23]. 

Another  possible  extension  is  to  model  the  service  channel  failures 
as  partial  breakdowns  due  to  sub-unit  failures.  Typically,  the  service 
channel  sub-units  can  be  removed  and  the  service  station's  ability  to 
perform  is  only  partially  affected.  An  example  of  this  is  an  electron- 
ics test  station  which  can  independently  repair  radio  and  radar  units. 
Failure  of  the  radio  section  may  not  affect  the  ability  to  service 
radars.  Thus  an  investment  in  service  channel  spares  should  be  con- 
sidered for  sub-components  and  not  for  entire  servers. 

This  situation  requires  a complex  reliability  analysis.  Spare  sub- 
units can  be  considered  as  elements  in  parallel  with  the  installed  sub- 
unit. Consequently,  the  entire  repair  station  could  be  modeled  as  a 
system  composed  of  sub-components  linked  in  series  and  parallel.  The 
purchase  of  an  additional  sub-unit  would  provide  a certain  marginal 
improvement  of  the  service  channel's  reliability.  If  this  phenomenon 
can  be  quantified,  then  performance  curves  could  be  derived  for  fixed 
levels  of  investment  in  server  sub-units  and  the  remaining  analysis 
should  be  similar  to  that  performed  in  Chapter  IV. 


Chapter  V suggests  one  final  area  for  future  research  with  its 
discussion  of  the  dearth  of  comparisons  between  the  different  time- 
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dependent  distributions  for  queues  found  in  the  literature.  It  would  be 
instructive  to  have  a standard  set  of  common  time-dependent  problems 
(including  non-stationary  inputs  to  the  queue)  to  provide  a means  for 
comparing  the  various  proposed  approximation  methods. 
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The  purpose  of  this  study  is  to  analyze  an  inventory,  maintenance 
system  for  recoverable  items,  that  is,  items  which  are  subject  to  repair 
when  they  fail.  The  repair  of  items  is  performed  by  a maintenance 
facility  which  has  a fixed  number  of  service  stations  or  channels  which 
are  also  subject  to  failure.  When  an  item  fails,  a demand  is  immediately 
placed  for  a like  replacement  from  a spare  pool.  The  failed  part  is  sent 
to  the  repair  facility  to  be  serviced  on  a first-come,  first-served  basis. 
The  spare  pool  is  replenished  when  repair  on  the  item  is  completed.  When 
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service  station  fails,  repair  is  initiated  immediately  and  the  failed  server 
is  replaced  by  an  operative  spare  server  if  one  is  available.  This  analysis 
is  limited  to  a single-echelon  system  with  no  outside  sources  of  supply  or 
repair. 


The  objective  of  this  study  is  to  model  the  system  described  in  order  to 
observe  the  relationship  of  system  performance  to  spare  stock  levels  and  service 
facility  design.  Specifically,  the  model  is  used  to  minimize  the  total  expected 
unit  backorders  given  an  investment  constraint  on  the  number  of  spare  items , 
service  channels  and  spare  servers  in  the  system.  For  long  range  planning 
purposes,  this  is  accomplished  for  a system  with  demands  which  are  stochastic 
and  stationary  in  nature.  An  extension  is  provided,  to  consider  the  case 
where  demands  are  non-stationary  and/or  the  time  dependent  behavior  of  the 
system  needs  to  be  described. 

expected  unit  backorders , a representation 
of  units  requiring  repair  is  needed. 
Approximations  ari  developed  us ing\dif fusion  techniques  (since  the  actual 
distributions  are/ difficult  to  express.  The  diffusion  approximation  is  applied 
to  an  optimization  problem  to  provide  the  best  allocation  of  investments  in 
the  system.  A simple  solution  algorithm  is  given. 

Finally,  a view  of  the  time-dependent  behavior  of  the  system  is  provided. 
The  -problem  is  decomposed  into  finding  the  distributions  for  (1)  the  number 
of  units  in  requiring  repair  given  no  service  channel  failures  and  (2)  the  time 
between  service  channel  failures.  We  provide  a brief  review  of  the  literature 
for  the  first  distribution  and  an  in-depth  study  of  the  latter  distribution. 


In  order  to  express  the  total 
for  the  distribution  of  the  number 


